Effect of Monovalent Ion Parameters on Molecular Dynamics

Article pubs.acs.org/JCTC

Effect of Monovalent Ion Parameters on Molecular Dynamics Simulations of G‑Quadruplexes Marek Havrila,†,‡ Petr Stadlbauer,†,§ Barira Islam,† Michal Otyepka,§ and Jiří Šponer*,†,‡ †

Institute of Biophysics, Academy of Sciences of the Czech Republic, Královopolská 135, 612 65 Brno, Czech Republic CEITEC - Central European Institute of Technology, Masaryk University, Campus Bohunice, Kamenice 5, 625 00 Brno, Czech Republic § Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacky University, 17. listopadu 12, 771 46 Olomouc, Czech Republic ‡

S Supporting Information *

ABSTRACT: G-quadruplexes (GQs) are key noncanonical DNA and RNA architectures stabilized by desolvated monovalent cations present in their central channels. We analyze extended atomistic molecular dynamics simulations (∼580 μs in total) of GQs with 11 monovalent cation parametrizations, assessing GQ overall structural stability, dynamics of internal cations, and distortions of the G-tetrad geometries. Majority of simulations were executed with the SPC/E water model; however, test simulations with TIP3P and OPC water models are also reported. The identity and parametrization of ions strongly affect behavior of a tetramolecular d[GGG]4 GQ, which is unstable with several ion parametrizations. The remaining studied RNA and DNA GQs are structurally stable, though the G-tetrad geometries are always deformed by bifurcated H-bonding in a parametrization-specific manner. Thus, basic 10-μs-scale simulations of fully folded GQs can be safely done with a number of cation parametrizations. However, there are parametrizationspecific differences and basic force-field errors affecting the quantitative description of ion-tetrad interactions, which may significantly affect studies of the ion-binding processes and description of the GQ folding landscape. Our d[GGG]4 simulations indirectly suggest that such studies will also be sensitive to the water models. During exchanges with bulk water, the Na+ ions move inside the GQs in a concerted manner, while larger relocations of the K+ ions are typically separated. We suggest that the Joung-Cheatham SPC/E K+ parameters represent a safe choice in simulation studies of GQs, though variation of ion parameters can be used for specific simulation goals.

■

relevancy, with K+ preferred over Na+.7,8,13,14 It is well established that Na+, due to its smaller size, can be observed in a variety of positions inside the GQ ionic channel. It can reside fully between two tetrads, but it can also be found closer to one of the tetrads, as well as in the center of one tetrad. We can even observe Na+ occupying all such positions in a single GQ structure.15,16 On the other hand, larger K+ is always placed equidistantly between two tetrads. Identity of the ions in the GQ channel may affect many of its properties, including the thermodynamic stability and in some cases even the dominantly populated topology.8,17 The broad spectrum of experimental methods employed in the investigation of GQs is complemented by theoretical approaches, mainly atomistic explicit-solvent molecular-dynamics (MD) simulations. For specific goals, quantum-chemistry (QM) calculations can be also applied in small model calculations.18−27 Contemporary MD is based on highly

INTRODUCTION Guanine quadruplexes (GQs) are specific secondary and tertiary structures, which can be formed by G-rich DNA and RNA sequences. GQs represent the most important noncanonical structure of DNA. Potential GQ-forming sequences are present at telomere ends1−3 and common in oncogenic promoters.4,5 There are many important regulatory functions in which GQs may be involved.6 The basic structural unit of the GQ is a guanine tetrad − a planar cyclic arrangement of four guanines stabilized by Hoogsteen base pairing (Figure 1). Tetrads stack together to form a 3D structure of GQ stems with the characteristic inner channel that binds cations (Figure 1). Monovalent ions play a fundamental role in formation and stabilization of GQs. The stabilizing effect can be ascribed to interaction of the cations with strong electronegative potential of guanine O6 oxygens oriented inward the channel.7−12 Besides the direct ion-tetrad coordination, the desolvation penalty of the ions is important for the overall thermodynamics balance. Among the wide range of ions stabilizing GQs, K+ and Na+ have the highest biological © 2017 American Chemical Society

Received: March 12, 2017 Published: June 28, 2017 3911

DOI: 10.1021/acs.jctc.7b00257 J. Chem. Theory Comput. 2017, 13, 3911−3926

Article

Journal of Chemical Theory and Computation

Figure 1. Examples of the simulated systems (starting structures) and a guanine tetrad. A) Full system of 1KF1 containing loops and the 5′ overhanging adenine. Residues inherent to tetrads and loops are blue and green, respectively. Cations are purple spheres, and O6 oxygens of guanines are red. B) Truncated d[GGG]4 system consisting only of guanine tetrads, i.e., with deletion of all loops and flanking nucleotides. C) Guanine tetrad with correct H-bonding pattern (dotted lines) and D) guanine tetrad with BHB. See Supporting Information Figure S1 for topologies of the other simulated systems.

approximate pair-additive nonpolarizable force fields (so-called molecular-mechanics description; abbreviated MM) which allow large-scale sampling of the conformational space. When using an appropriate DNA force field, GQs stabilized by monovalent ions are structurally stable in MD simulations.28,29 This indicates that the overall electrostatic stabilization of the GQs by monovalent ions is very well captured by the MD technique; for a detailed overview see Sponer et al.28 The ionGQ electrostatic stabilization is indeed widely considered as the dominant energy term stabilizing GQ molecules.8,30,31 GQs are among the best-behaving systems in nucleic acids simulations, which can be related to the structural stiffness of the GQ stems.29 Despite this, already the first MD simulations indicated occurrence of more subtle inaccuracies in the description of the ion-GQ interactions, stemming from the simplicity of the force field.30 Some of the problems may be related to specific ion parametrizations, but most of them are due to the general approximation of the MM description. Contemporary force fields are based on constant point charge approximation, which neglects the electronic structure and consequently polarization effects. This approximation renders ions as spheres differing in van der Waals (vdW) parameters, mass, and point charge. Therefore, interaction potentials of monovalent ions differ only in their vdW parameters, which is not sufficient to tune the ion parameters simultaneously for different interactions and cation-binding patterns. The following specific problems with quantitative description of the cations inside the GQ stems have been identified. First, sodium cations during MD are typically located out of the plane of tetrads, in contrast to the X-ray structures.8,15 This may be partially due to specific ion parametrizations, as the ion parameters are typically fitted to reproduce solvation energies in water and are thus not optimized for interactions with Gtetrads. Another artifact is the appearance of bifurcated hydrogen bonding (BHB; see Figure 1), somewhat deformed geometries of the tetrads which have never been observed in atomistic experimental structures.30 Although the BHB of tetrads does not cause any substantial perturbations of the simulated structures, it indicates that the ion-tetrad interactions are not fully balanced. Probably the most serious problem in

the description of ions inside GQs has been visualized recently using benchmark QM calculations for interactions between multiple cations inside the GQ stems.21 When two or more ions are present inside the stem, their effective mutual repulsion is significantly overestimated by MM. Real ions electronically communicate with each other through the stacked tetrads via a complex relay of polarization effects. Thus, the mutual electrostatic repulsion between the cations is typically softened compared to interaction of two +1 point charges. In reality, any change of the coordinates of the ions and tetrads modulates the electronic structure. Such continuous change of electronic structure cannot be captured by fixed-charge models. This limitation of force fields is likely responsible for inability of MD simulations to keep bound cations at the stem-loop junctions of the d[G4T4G4]2 diagonal loop GQ.32 These ions are expelled due to their excessive repulsion with the adjacent channel ions.19,21 Stabilization of junction cations was reported when using recalculation of the DNA charges specifically for the folded d[G4T4G4]2 GQ conformation.19 However, such an adapted force field still remained pair-additive with fixed charges and was thus specific for one structure.19 Perhaps, the description of the interactions between GQs and their bound cations could be in the future improved using polarizable force fields, which are under intense development.33−35 However, their tuning is difficult, and they are not yet matured for routine nucleic acids simulations. Thus, it is very likely that atomistic simulations with the current class of pair-additive force fields will dominate computational studies of GQs in the foreseeable future.36−54 We therefore need to understand their limitations and minimize errors. Despite that it is generally proven that with good force field the dominant electrostatic stabilization of GQs by monovalent ions is well captured, there has been no systematic assessment of the available cation parametrizations for long timescale simulations of GQs. Several cation parametrizations have been commonly used in the past, and it is not clear if their performance in GQ simulations is identical or not. Because the ion parameters are not parametrized directly for simulations of GQs, it is not a priori guaranteed that those parameters that are generally recommended for biomolecular simulations are also 3912


Article

Journal of Chemical Theory and Computation the best ones for the description of the ions inside the GQs. Here we provide extensive testing of different parametrizations of NaCl and KCl salts in GQ simulations and address several questions associated with the effect of ion parametrization on properties of GQ in MD simulations. (i) How do the ion identity and parametrization affect structural stability of GQs in MD simulations? (ii) What is the mobility of different ion types in the channel and frequency of their exchange with the bulk? (iii) How common is the BHB geometry in MD simulations? (iv) Which of the tested cation parameters is the most appropriate for MD simulations of GQs? Regarding the last point, we wish to underline that due to the basic force-field approximation noted above, we cannot obtain an error-free description of the cation dynamics and interactions inside the GQ stems. However, our goal is to select parameters which minimize inaccuracies in the simulations. The present study demonstrates large differences of the available monovalent ion parametrizations on the cation−tetrad interactions. Simultaneously, we show that there is no detectable effect of the tuning of DNA backbone dihedral potentials on the ion-tetrad interactions and ion dynamics. The present manuscript is based on analysis of 580 μs of explicit-solvent MD simulations, which in our opinion is the largest total accumulated simulation time in a single GQ simulation study. The individual simulations are either 5 or 10 μs long.

Table 1. Lennard-Jones Parameters of Tested Monovalent Ions ion modela Na+ JC AA IOD HFE K+ JC DK IOD HFE Cl− JC SD IOD HFE

ε (kcal/mol)

rm/2 (Å)

0.3526 0.0028 0.0291 0.0264

1.212 1.868 1.465 1.454

0.4297 0.1000 0.1702 0.1513

1.593 1.870 1.745 1.683

0.0128 0.1000 0.5315 0.6437

2.711 2.469 2.162 2.308

Parameters of ions are abbreviated as follows: AA − AMBER-adapted Aqvist, JC − Joung and Cheatham, DK − Dang and Kollman, SD − Smith and Dang, HFE − Li, Song, and Merz ions optimized for hydration free energy reproduction; IOD − Li, Song, and Merz ions optimized for ion−oxygen distance reproduction. Note that SPC/E related parameters are shown. a

■

METHODS Description of Ions in Force Fields. In most common force fields, the interactions of ions are described by a simple intermolecular potential, which consists of Coulombic and Lennard-Jones terms, representing electrostatic and vdW interactions, respectively. Use of the Coulombic potential is facilitated by approximating the electronic structure of molecules by fixed point charges located at the atom nuclei. The Lennard-Jones potential υLJ depends on the distance between interacting particles r and two parameters: ε − the well depth of potential, and rm − the distance at which the energy is a minimum. rm can be substituted by the collision diameter σ (σ = rm/21/6), the distance for which energy is zero. Lorentz− Berthelot combining rules are employed in the AMBER software package to obtain effective vdW parameters for pairs of particles.

Plots of ion−ion and ion−oxygen LJ potentials used in this study are in Figure 2. Note that some additional ion parameters which were used only in a limited series of calculations will be specified in the respective part of the Results section. Simulated Systems. We simulated a variety of DNA GQs based on the human telomeric sequence d(TTAGGG)n, the human telomeric RNA sequence r(uuaggg)n (TERRA) GQ, and the c-kit promoter DNA GQ (Table 2). Schematic representation of all simulated systems is available in Supporting Information Figure S1. We divided the simulations into three groups based on the primary purpose of the simulations (Table 3). Group I is based on monomolecular DNA parallel GQ with three tetrads (d[AGGG(TTAGGG)3], PDB code 1KF161). This system has been chosen because it contains only propeller loops, which are located on the sides of the GQ core and cannot close the central channel by stacking onto the peripheral tetrads as is often seen in various other folds (see Figure 1A). An advantage of this system is also straightforward comparison to RNA, which forms exclusively parallel GQs. The Group I set of simulations has been designed for general study of the influence of ion parametrization on GQ. It contains simulations with all tested salt parametrizations and their combinations. Besides the full 1KF1 system, we prepared a truncated 1KF1 system without loops, i.e., d[GGG]4 (see Figure 1). Simulations of the native system are designated with “F” (Full) in the simulation name, while those of the truncated system contain “T” (Truncated) in the name (Table 3). The d[GGG]4 system is very suitable for our study, since it is more dynamic than the complete system, and the effects of different ion parametrizations are better manifested. Each system under given conditions was simulated in three independent MD runs to improve sampling and to prevent biasing of results by random rare events. Group II is based on bimolecular RNA parallel GQ with three tetrads and two propeller loops (r[uaggguuagggu]2, PDB code 3IBK62) and is used to characterize RNA GQs. Only a

⎡⎛ σ ⎞12 ⎛ σ ⎞6 ⎤ υLJ(r ) = 4ε⎢⎜ ⎟ − ⎜ ⎟ ⎥ ⎝r⎠ ⎦ ⎣⎝ r ⎠

In this study we tested monovalent ion parameters that are widely used in the MD community for nucleic acids simulations (see Table 1): Aqvist55 (AA) Na+ with Smith and Dang56 (SD) Cl−, which were default NaCl parameters in older AMBER versions; Dang and Kollman57 (DK) K+ with SD Cl−, a widely used combination for KCl in older simulations; Joung and Cheatham58 (JC) NaCl and KCl, improved parameters popular in recent studies; and recently developed Li-Song-Merz59 HFE (optimized to reproduce the experimental hydration free energies) and IOD (optimized to reproduce experimental ion−oxygen distance) ion sets. Ion parameters optimized for the SPC/E water model were used (see below). Since mass and charge of certain monovalent ion are fixed, the parametrizations differ only in ε and rm values. Overview of existing parameters can be found in a paper by Joung and Cheatham58 or in the reference manual to the AMBER software package.60 A short comparison can also be found in Li et al.59 3913


Article


Group III contains a variety of DNA GQs with different folds based on structures denoted by PDB codes 1K8P,63 143D,64 2KF8,65 2KM3,66 2HY9,67 2JPZ,68 2MBJ,69 and 3QXR70 (see Table 2 and Supporting Information for details). Group III simulations are done with three different DNA force fields but only with a limited number of ion parametrizations. The aim of these simulations is to verify whether the data obtained in Group I and II simulations are generally valid and to examine if details of the nucleic acid force-field backbone parametrization may affect the results. Simulation Protocol. Explicit-solvent MD simulations were performed using the AMBER14 software package.60 The xLEaP and parmed.py modules of AMBER14 were used to prepare starting topologies and coordinate files. The following variants of the AMBER nucleic acid force field were used: (i) bsc0χOL371−73 for all RNA simulations, (ii) OL15 force field,50 which consists of parm9971 modified by bsc0,72 χOL4,37 εζOL1,74 and βOL150 corrections for Group I simulations; (iii) bsc0χOL4, bsc0χOL4 + εζOL1, and OL15 for Group III simulations. For an overview of the used force fields see Sponer et al.28 and Galindo-Murillo et al.75 OL15 represents a complete reparametrization of the DNA backbone dihedrals compared to the 1995 Cornell et al. AMBER force field71 and is at the performance limits that can be achieved by tuning of the DNA backbone dihedral potentials. Note that all the tested versions use the original Cornell et al. nonbonded terms, and thus any effect of the dihedral force-field variant on ion-tetrad interactions could be only indirect. Starting structures were immersed in a truncated octahedral box of SPC/E explicit water molecules,76 and ions were added to obtain 0.2 M excess salt concentration. Most of the ion models are parametrized for specific water models. Therefore, we used ion parameters optimized for the SPC/E water model. A small set of additional simulations was done with either OPC77 or TIP3P78 water models, which are discussed separately at the very end of the Results and Discussion section. Thus, unless specified otherwise, all presented data are with the SPC/E model and corresponding ion parameters. Prepared systems were minimized, heated, and equilibrated in several steps before each MD run. The velocities of atoms were randomized in the first equilibration run to get independent starting structures for simulations with otherwise identical properties. Hydrogen mass repartitioning was applied together with the SHAKE algorithm79 to allow for the 4 fs integration time step.80 Simulations were performed using the particle mesh Ewald method.81,82 We used the standard

Figure 2. Different LJ potentials used in this study. Top − ion−ion interaction, bottom − ion−O6 interaction.

subset of ion parametrizations has been investigated. The native 3IBK system contains overhanging nucleotides, which are involved in interactions with another GQ in the crystal cell. Since these overhangs lacked an interacting partner and created undesired contacts in our preliminary simulations, we removed them, and our “full” system refers to a system without the overhangs, i.e., r[ggguuaggg]2. We also simulated its truncated variant r[ggg]4. Two independent MD runs were performed for each system under given conditions. Table 2. Simulated Systems name

sequence

no. of chains

no. of tetrads

strand orientation

loops

PDB

T_1KF1 F_1KF1 T_3IBK F_3IBK 1K8P 143D 2KF8 2KM3 2HY9 2JPZ 2MBJ 3QXR

[d(G3)]4 d[A(G3T2A)3G3] [r(G3)]4 [r(G3U2AG3)]2 [d(TAG3T2AG3T)]2 d[AG3(T2AG3)3] d[G3(T2AG3)3T] d[AGGG(CTAGGG)3] d[A3G3(T2AG3)3A2] d[(T2AG3)4T] d[(T2AG3)4T2A] d[AG3AG3CGCTG3AG2AG3]

tetramolecular monomolecular tetramolecular bimolecular bimolecular monomolecular monomolecular monomolecular monomolecular monomolecular monomolecular monomolecular

3 3 3 3 3 3 2 3 3 3 3 3

parallel parallel parallel parallel parallel antiparallel antiparallel antiparallel (3+1)hybrid (3+1)hybrid antiparallel parallel

none 3x propeller none 2x propeller 2x propeller 1x diagonal; 2x lateral 1x diagonal; 2x lateral 3x lateral 1x propeller; 2x lateral 1x propeller; 2x lateral 1x propeller; 2x lateral 3x propeller; 1x lateral

1KF1 1KF1 3IBK 3IBK 1K8P 143D 2KF8 2KM3 2HY9 2JPZ 2MBJ 3QXR

3914


Article

Journal of Chemical Theory and Computation Table 3. List of Simulations simulation namea

run length (ns)

no. of equivalent runs

5000 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000

3 3 3 3 3 3 3 3 3 3 2

5000 5000 5000 5000

Group I 1 T_1KF1_JC_NaCl 2 T_1KF1_JC_KCl 3 T_1KF1_AA_NaCl 4 T_1KF1_DK_KCl 5 T_1KF1_IOD_NaCl 6 T_1KF1_IOD_KCl 7 T_1KF1_HFE_NaCl 8 T_1KF1_HFE_KCl 9f T_1KF1_JC_Na/KCl 10f T_1KF1_JC_K/NaCl 11f T_1KF1_JC_AA_K/NaCl 12 F_1KF1_JC_NaCl 13 F_1KF1_JC_KCl 14 F_1KF1_AA_NaCl 15 F_1KF1_DK_KCl Group II 16 T_3IBK_JC_NaCl 17 T_3IBK_JC_KCl 18 T_3IBK_AA_NaCl 19 F_3IBK_JC_NaCl 20 F_3IBK_JC_KCl 21 F_3IBK_AA_NaCl Group III 22−48e 143D; 1K8P; 2KM3; 2HY9; 2JPZ; 2KF8; 2MBJ; 3QXR; F_1KF1

conc. (M)c

ionsb

force fieldd

0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2

OL15 OL15 OL15 OL15 OL15 OL15 OL15 OL15 OL15 OL15 OL15

3 3 3 3

NaCl (JC) KCl (JC) NaCl (AA) KCl (DK) NaCl (IOD) KCl (IOD) NaCl (HFE) KCl (HFE) Na+/ KCl (JC) K+/NaCl (JC) K+(JC)/NaCl (AA) NaCl (JC) KCl (JC) NaCl (AA) KCl (DK)

0.2 0.2 0.2 0.2

OL15 OL15 OL15 OL15

5000 5000 5000 5000 5000 5000

2 2 2 2 2 2

NaCl (JC) KCl (JC) NaCl (AA) NaCl (JC) KCl (JC) NaCl (AA)

0.2 0.2 0.2 0.2 0.2 0.2

bsc0χOL3 bsc0χOL3 bsc0χOL3 bsc0χOL3 bsc0χOL3 bsc0χOL3

10000 10000 5000

1 1 1

NaCl (JC) NaCl (JC) KCl (JC)

0.2 0.2 0.2

bsc0χOL4 bsc0χOL4 + εζOL1 OL15

The letter T in the name of the simulation stands for “Truncated”, while the letter F stands for “Full” system. bAbbreviations of ion parameters are in parentheses. AA − AMBER adapted Aqvist Na+ with Smith and Dang Cl−; DK − Dang and Kollman K+ with Smith and Dang Cl−. The rest of the abbreviations follows Table 1. All simulations specified in this table were done with the SPC/E water model. cThe listed value represents total concentration of anions (i.e., excess salt). dOL15 stands for parm99 modified by bsc0, χOL4, εζOL1, and βOL1 corrections. eFor each of the nine listed systems, we performed three simulations with different parameters, as indicated in the three rows (except for the 2KM3 system, which was simulated only in OL15 and bsc0χOL4). The F_1KF1 simulation in OL15 is additional to the simulations of row 13. fDifferent parameters were used for ions in the channel and in bulk. Channel ions are at the left side of the slash, while bulk ions are at the right side. a

simulations were somewhat less stable. Occasional reversible vertical strand slippages were observed (Figure 4); however, the prevailing structure was a stable GQ. The vertical strand slippage drives a GQ to a structure containing two G-tetrads, a G-triad and an overhanging G. Bases in the triad can flip from anti to syn conformation (see below), which then stabilizes the slipped structure, as the syn nucleotide prevents the backslippage (see Supporting Information Figure S2).38 Such stabilization of the slipped structure was observed in one of three performed simulations. Strand-slippage is a mechanism which has been suggested as a key component of formation of parallel-stranded monomolecular GQs,48 and this type of movement may be also relevant for folding and unfolding of parallel-stranded tetramolecular GQs.38 To the best of our knowledge, it is for the first time that a spontaneous strand slippage has been achieved in simulations of an ion-occupied GQ. This reflects the length of our simulations as well as some structural volatility of the d[GGG]4 construct. Different disruption of the GQ was observed in T_1KF1_IOD_KCl simulations. One of the three simulations developed a major disruption of the 5′-terminal tetrad, which was followed by vertical strand slippage leading to a structure containing two G-tetrads, a G-triad and an overhanging G. Gtriad then flipped to syn conformation and stabilized the slipped

simulation setup that can be found in ref 83, and its details are given in the Supporting Information.

■

RESULTS AND DISCUSSION

Ion Type and Parametrization Dramatically Affect Structural Stability of the d[GGG]4 System, while Stability of the Full 1KF1 System Is Affected Negligibly. We observed large differences in structural stability of the d[GGG]4 system under various ionic conditions. By stability of the simulated system in the context of our study we mean its overall structural behavior (with respect to the starting structure) in the simulations, including fluctuations reflected by RMSd and mainly occurrence of larger structural changes and perturbations that are described below. Note that the presently studied systems are sufficiently simple that we can obtain a good estimate of how a given system behaves from its basic description. Representative RMSD plots are shown in Figure 3. All T_1KF1_JC_KCl simulations were perfectly stable with average RMSD to the starting structure around 1.3 Å. No fluctuations or disruptions of the system were observed. T_1KF1_HFE_KCl simulations were also relatively stable, and only one simulation in this set showed a single event of reversible disruption of the 3′ tetrad. T_1KF1_JC_NaCl 3915


Article


Figure 3. Development of selected properties of GQ in MD simulations of the truncated 1KF1 d[GGG]4 systems. One representative trajectory is displayed for each condition. KCl simulations are in the left column, NaCl simulations are in the middle column, and simulations which have different ions in the channel and in bulk are in the right column. Each graph consists of two parts. The upper part captures development of RMSD of atom positions (gray) and diagonal distances of O6 oxygens (in black and green; the inner tetrad is visualized, see below for an explanation of the significance of the distances). The bar graphs in the lower part monitor exchanges of ions in the channel. a and b represent two typical positions of ions in the channel of GQ. a stands for the position between the 5′ tetrad and the inner tetrad, while b represents the position between the inner and 3′ tetrad. Distinct ions are represented by different colors, and any change in the color indicates an exchange of an ion in the channel. Note that the description of ions in the channel is reliable only in stable systems with RMSD lower than 2 Å. Our algorithm analyzes geometrical parameters even if one tetrad is lost or if the GQ is disrupted, but the results of such analyses are already less relevant. For a complete dataset see Supporting Information Figure S3.

structure (similarly to the mentioned T_1KF1_JC_NaCl simulation). Another simulation developed a small reversible disruption of the 3′ tetrad interrelated with a temporary loss of one channel ion. In comparison to T_1KF1_JC_NaCl simulations where we primarily observed strand slippages, the disruption of the outer tetrad is primary perturbation in the T_1KF1_IOD_KCl simulations. The third simulation did not develop any disruptions. The preference of 5′ terminal guanines for the syn orientation has been documented in several studies.20,29,84,85 The 5′ terminal guanine is unique by the presence of the free O5′-H terminus, which can form intramolecular O5′H···N3 Hbond in syn conformation. This hydrogen bond affects the equilibrium of syn versus anti conformation both in computations and in experiments, by forcing the 5′-terminal

Figure 4. Schematic representation and the structure of the d[GGG]4, in which one strand underwent the vertical strand slippage.38,48 The slipped strand is green, the rest of the guanines are blue, and guanine O6 oxygens are red.

3916


Article


Figure 5. Development of studied properties in MD simulations of full 1KF1 systems. One representative trajectory is displayed for each condition. For more details see the caption of Figure 3. For a complete dataset see Supporting Information Figure S4.

Figure 6. Development of selected properties in MD simulations of RNA GQ based on the 3IBK system. One representative trajectory is displayed for each condition. Truncated r[ggg]4 and full systems are in the upper and lower rows, respectively. For more details see the caption of Figure 3. See Supporting Information Figures S5 and S6 for the full dataset.

stable until 4.3 μs, but considerable structural deterioration finally appeared at the end of the simulation, which likely would lead to a complete loss of GQ upon prolongation. Even faster loss of the GQ structure and formation of random aggregates were observed in all three T_1KF1_HFE_NaCl simulations which can be considered as the least stable. In summary, the differences in ion parameters had unexpectedly large impacts on the simulation behavior. Since the d[GGG]4 has not been experimentally characterized (due to technical difficulties, such as multimerization) and its thermodynamic stability is not known, we cannot tell which of the observed simulation behaviors is more correct. The results nevertheless show the clear influence of ion parameters on the studied system. Considering the generally known stability of GQs, we hypothesize that those ion parametrizations that show profound instability of d[GGG]4 on the present simulation timescale are likely less suitable for GQ simulations. It should be noted that our set of three simulations for each ion parameter set was not sufficient to provide a quantitative statistics of large structural perturbations, which in addition was not our main goal. The performed simulations should be nevertheless capable of qualitatively depicting behavior of the d[GGG]4 simulated with different cations in the channel. To confirm that the observed differences in d[GGG]4 simulations were caused by ions in the GQ channel and not by the bulk ions, we ran simulations with different parametrization of ions in the channel and in bulk. Simulations with

guanines into syn conformation.29 It has been first predicted20,29 and then observed in NMR86 and X-ray87 structures. The effect of O5′-H terminus is maximized in our truncated system by four 5′ terminal guanines, and it leads to the swift stabilization of the syn conformation of the 5′ triad in the above-noted T_1KF1_IOD_KCl simulation. Note that the transition to syn from the starting all anti conformation requires that the 5′ tetrad is first disrupted, allowing for anti-syn transition of guanines. Both strand slippage and disruption of the 5′ tetrad were observed in T_1KF1_DK_KCl simulations. Disrupted structures were prevailing; however, at least two complete tetrads were always present, and the GQ was not fully lost (Supporting Information Figure S2). Alteration of glycosidic orientation of 5′-terminal guanines was often present. T_1KF1_IOD_NaCl simulations developed large disruptions of the GQ core. In two of the three simulations strand slippages initiated structural changes that developed to almost unstructured aggregates of DNA strands centered around a trapped Na+ ion, which were however able to refold back to systems with two tetrads. The prevailing structure was a somewhat disrupted GQ with two complete tetrads. The third simulation led to random aggregates and stacked structures. T_1KF1_AA_NaCl simulations showed very low stability. Complete loss of the GQ and presence of random aggregates and stacked structures were present in two simulations (see Supporting Information Figure S2). The third simulation was 3917


Article


Figure 7. Fluctuations of backbone torsions in truncated DNA and RNA GQs. The graphs represent time development of all backbone torsions in the two consecutive backbone suites89 of one chain of the truncated systems simulated in JC KCl. Larger fluctuations of the DNA backbone (d[GGG]4 − upper graphs) in comparison with RNA (r[ggg]4 3IBK − bottom graphs) are visible.

Figure 8. Mechanisms of ion exchanges observed in simulations of GQ. The ions have different colors to better visualize the exchanges: A) concerted movement of ions in the channel observed for various Na+ parametrizations; B) blockage of ion exchange by the overhanging adenine (depicted as black rhomboid) interacting with O6 oxygens of the 5′ tetrad in the full 1KF1 system; C) loss of Na+ and possible insertion of a water molecule in T_1KF1_IOD_NaCl simulations; D) ion exchanges in simulations with IOD and DK K+ including presence of water molecule.

K+ in the channel and Na+ in bulk were performed in two combinations. In the T_1KF1_JC_K/NaCl system, JC parameters were used for all ions, while in T_1KF1_JC_AA_K_NaCl, JC parameters for K+ and AA parameters for Na+ were used. The third combination, designated as T_1KF1_JC_Na/KCl, consisted of Na+ in the channel and K+ in bulk employing JC parameters. This cross-comparison verified that the observed effects were linked to the channel

ions. The best verification was high stability of the system with JC K+ in the channel and AA Na+ in bulk, because the system simulated solely with JC KCl was stable, while the system with AA NaCl was very unstable. Moreover, the T_1KF1_JC_K/ NaCl simulations were similar to the T_1KF1_JC_KCl ones, and the T_1KF1_JC_Na/KCl simulations were similar to the T_1KF1_JC_NaCl ones, until Na+ left the channel and K+ was 3918


Article


balance of cation−tetrad, cation−water, and cation−cation interactions. The simulations T_1KF1_JC_KCl, in which the GQ was perfectly stable, did not show any exchanges of ions between the channel and bulk (Figure 3). However, K+ ions could jump out of the channel, stay centered close to the outer tetrad, and return to its original position inside the GQ. The whole process occurred on the timescale of several ns and appeared a few times in a 5 μs long simulation. Similarly, the T_1KF1_HFE_KCl simulations did not show any exchanges except of a short time window with a disrupted 3′ tetrad in one simulation. In striking contrast, the simulations T_1KF1_JC_NaCl, which also largely preserved intact GQ, exhibited rapid ion exchange with residency time of the ion in the channel from tens to hundreds of ns. Thus, structural stability of d[GGG]4 GQ did not correlate with the mobility of ions in the channel. In addition, the Na + ions in the channel of the T_1KF1_JC_NaCl system showed correlated movements. More accurately, they were moving as a couple (or as a triple, when a third ion approached from the bulk) of ions. When one ion left its position, another ion replaced it simultaneously. The exchange process always started by approach of a cation from the bulk to proximity of an outer tetrad (either 5′ or 3′), where the ion was caught nearby the GQ channel. Afterward, the three ions (two in the channel and one in its proximity) moved simultaneously through the channel so that the outer ion entered the channel and the farther inner ion left the channel at the opposite side (Figure 8A), so it looks like the incoming cation was “catalyzing” the ion exchange event. This coordinated movement explains why the presence of the 5′ overhanging adenine in 1KF1 interacting with O6 oxygens of the 5′ tetrad impaired ion exchange despite the fact that the channel remained open at the 3′ end (Figure 8B). Similar behavior of Na+ ions was observed also in T_1KF1_IOD_NaCl and T_1KF1_AA_NaCl simulations; however, the exchanges were much slower. Movement of ions in the channel was correlated, and the mechanism of ion exchange described above prevailed. However, in the T_1KF1_IOD_NaCl simulation, we could sometimes observe simultaneous movement of two ions in the channel even without the presence of a third cation in its proximity, leading from time to time to a GQ with one ion in the channel and one in its proximity. The outer ion could then be lost, temporarily resulting in a GQ with a single ion in the channel (see Figure 8C). A water molecule can enter the channel and occupy the empty position in such a case. The T_1KF1_IOD_KCl and T_1KF1_DK_KCl simulations showed a different mechanism of ion exchange with a lower exchange rate. The K+ ions moved independently. One ion was able to leave its position in the channel, while the other one remained at the same place. As a result, a GQ with a single ion in the channel was present until another ion entered it (see Figure 8D). The process of binding a new ion into the empty position took from a few to hundreds of ns. Two factors that affect the speed of this process were the rate at which bulk ions sampled the position close to the channel and entry of a water molecule to the empty position in the GQ channel, which needed to be expelled first. The fact that the Na+ ions exchange with the bulk differently than the K+ ions is not surprising. The larger K+ ion faces a sizable barrier to cross the tetrad plane, which is of a steric (short-range repulsion) origin. This naturally reduces the

soaked in. This is an entirely expected result due to the dominant role of the channel cations in GQ−ion interactions. We then ran the simulations of the full 1KF1 GQ in several combinations of ion types and parameters (see Table 3). In sharp contrast to the simulations of the truncated system, simulated full parallel GQs with loops were stable, and we did not see any considerable differences even between JC KCl and AA NaCl ion conditions, which showed an extremely different behavior for the d[GGG]4 system (see Figure 5). Stability of the full GQ was similar under all the conditions, revealing a major stabilizing role of the loops, which were likely attenuating the fluctuations of the whole system and thus reducing the likelihood of the ion exchange, strand-slippage, and other structural perturbation. Parallel RNA G-Quadruplex Is More Structurally Stable than Its DNA Equivalent. The RNA GQ was simulated in JC KCl, JC NaCl, and AA NaCl conditions. Both full and truncated systems were stable in all tested salts, which was strikingly different from the unstable truncated DNA GQ (Figure 6). The effect of the RNA’s 2′OH group on structural stability of the GQs is 2-fold. The OH group forms hydrogen bonds with acceptors within the RNA GQ molecule, directly increasing its stability.62,88 The identity of the acceptor atom depends on the exact sugar pucker. The presence of the 2′OH group leads to a preference for the C3′ endo sugar pucker, which is native for the RNA parallel GQ. Different sugar pucker in 3IBK in comparison with 1KF1, which has exclusively C2′ endo sugar pucker, is compensated by correlated changes of backbone torsion angles, allowing for the same orientation of nucleobases in both GQs. The origin of the larger structural stability of RNA GQ (compared to DNA) in our simulations can be attributed to both mentioned effects. While the RNA backbone experienced very small fluctuations of torsion angles around their mean values (Figure 7), larger fluctuations of torsion angles were observed in the DNA simulations (Figure 7). These fluctuations of the DNA backbone were present not only in the unstable d[GGG]4 simulations but also in the stable T_1KF1_JC_KCl system, pointing to intrinsic flexibility of the DNA backbone in simulations of parallel DNA GQ. Mobility of Ions in the GQ Channel Depends on Ion Type and Parameters and Does Not Correlate with the Overall GQ Structural Stability. We have evaluated exchanges of ions between the channel and bulk for all simulated systems. We analyzed those parts of simulations, where GQs with intact tetrads were present. Moreover, only Group I of simulations (see Methods) contained systems with all tested salt conditions, so the following considerations are based mainly on the d[GGG]4 system. Further, in the full 1KF1 system, the 5′ overhanging adenine often formed H-bonds with O6 oxygens of the 5′ tetrad and blocked the approach of ions to the 5′ side of the channel from the bulk, resulting in the slowdown of ions exchanges (see Figure 8). Order of ion mobility in the d[GGG]4 GQ in different ionic conditions was approximately KCl JC < KCl HFE < KCl IOD < NaCl AA < KCl DK < NaCl IOD < NaCl JC (Figure 3; Supporting Information Figure S3). The three simulations with NaCl HFE ions did not provide sufficient fraction of stable parts of trajectories suitable for analysis of ion exchanges. Moreover, we noticed that different ion parametrizations resulted in systematically different types of exchange events, revealing visible impact of the specific parametrizations on the 3919


Article


Figure 9. X-ray and bifurcated geometries of tetrads and an example graph capturing different tetrad geometries: A) X-ray geometry of guanine tetrad; B) symmetrically bifurcated guanine tetrad; C) asymmetrically bifurcated guanine tetrad; D) graph explaining the representation of different inner tetrad geometries in the figures in this paper. Diagonal distances of O6 oxygens are in green and black, respectively. RMSD of all GQ heavy atoms is in gray. Region A + B depicts fluctuations of the system between X-ray and bifurcated geometry, and region B represents symmetrical bifurcation of the system. There are two possible conformations of tetrad with asymmetric bifurcation resulting from symmetry of the tetrad. Region C represents asymmetrical bifurcation locked in one conformation, while C* represents dynamical asymmetrical bifurcation with the two possible conformations quickly alternating.

3.5 Å, the distance between two JC K+ ions was ∼3.65 Å. Na+ was more dynamic and was able to visit a plane of tetrad, but the population of the in-plain Na+ was less than 1%. On the contrary, dynamics of K+ in the channel was visibly restricted, and the ions were literally sitting at their positions. Bifurcated Hydrogen Bonding of Tetrads Is Omnipresent and Strongly Affected by Ion Types and Parameters. BHB of tetrads is an artifact occurring in MD simulations of GQs, which was first described by Spackova et al.30 as a local difference in cyclic hydrogen bond geometry of tetrads. In the originally reported BHB variant, guanines of inner tetrads were shifted in such a way that BHBs were formed between the N1 nitrogen of one guanine and the N7 and O6 atoms of an adjacent guanine. This artifact was originally observed in telomeric DNA GQ of Oxytricha nova d[G4T4G4]2 and in the tetrameric d[GGGG]4 quadruplex simulated in the TIP3P water model and Amber-adapted Aqvist ion parameters, specifically for the inner tetrads.30 However, it is present to a variable extent in a variety of GQ topologies simulated with diverse water and ion models.19,88,90 In the present simulations, we also observed alternative variants of the BHB. While BHB described in the literature was symmetric (see Figure 9B) with O6 oxygens being corners of a virtual square, we also observed asymmetric geometry of tetrads with BHB resulting in rhombic geometry of O6 oxygens (see Figure 9C). The change of geometry of O6 oxygens related to BHB is reflected by the distances of two opposite O6 oxygens. These distances are lower in comparison with X-ray values in a symmetric BHB tetrad. The asymmetric BHB is characterized by one diagonal O6−O6 distance much shorter and the other one longer in comparison to the X-ray structure (see Figure 9). To analyze the presence of BHB, we evaluated only those parts of simulations, where stable GQs containing intact tetrads were present. In accordance with Spackova et al.,30 we observed BHB mainly in the inner tetrad. The scale of described perturbations was affected by ion type and parametrization; however, BHB was present to a certain extent in all our simulations. In general, the simulations with KCl showed symmetric BHB, while in all the NaCl simulations the asymmetric BHB geometry was present. JC KCl ions showed a very consistent effect on all simulated 1KF1 and 3IBK systems. Outer tetrads preserved X-ray conformation, and the inner tetrad exhibited symmetric BHB,

probability that two K + ions cross the tetrad planes simultaneously, in contrast to the correlated movements of the Na+ ions. This feature is likely correctly reflected by the force fields, since it is in agreement with benchmark QM calculations,21 albeit the later calculations could have been done only as one-dimensional potential energy scans not taking into account thermal fluctuations of the tetrads. However, besides the difference between the Na+ and K+ ions, the present simulations reveal a surprisingly large sensitivity of the ion behavior to the specific van der Waals parametrizations. Simulations of RNA GQ and Full DNA GQ Confirm the Trends of Ion Movements Observed for the Truncated DNA GQ. The simulations of the full 1KF1 system and the truncated 3IBK RNA system showed a very low number of ion exchanges (Figures 5 and 6). It can be ascribed to the presence of the 5′ overhanging adenine in the full 1KF1 system (see above) and higher overall stability (reduced structural fluctuations) of both systems in comparison to d[GGG]4. The increased stiffness likely reduced the likelihood of crossing the energy barriers for larger relocations of the ions. Although fewer simulations were done on these systems rendering a number of observed ion exchanges insufficient for quantitative analysis, the basic assessment of ion dynamics in these systems should be reliable. All ion exchanges observed in the full and truncated 3IBK and the full 1KF1 systems proceeded similarly to the d[GGG]4 simulations with the corresponding ion types. Na+ displayed the correlated movement, while K+ showed separation of the cation movements leading to the formation of a GQ with a single ion in the channel. No exchanges were observed in F_1KF1_JC_KCl and T_3IBK_JC_KCl systems, while one exchange occurred in F_3IBK_JC_KCl, T_3IBK_JC_NaCl, T_3IBK_AA_NaCl, and F_1KF1_AA_NaCl systems. Several exchanges were observed in F_1KF1_DK_KCl, F_1KF1_JC_NaCl, F_3IBK_AA_NaCl, and F_3IBK_JC_NaCl simulations. This observation allows us to conclude that the described ion behavior is valid for simulations of all parallel GQs in general. Movement of smaller Na+ ions in the channel of GQs is concerted, while bigger K+ prefers to move independently. The positions of Na+ and K+ in the GQ channel during the course of the simulations were not equivalent. While the average separation of two JC Na+ ions in the channel was 3.4 to 3920


Article


the ion dynamics inside the GQ stem. However, we consider this situation as quite unlikely and suggest that ion dynamics analyzed in this work are not affected by the backbone dihedral potentials. In other words, we assume (and the data described in the subsequent paragraph are consistent with this assumption) that the results of our work are valid for all DNA force-field versions irrespective of the backbone dihedral potentials, provided they guarantee a stable description of the GQ stem. This, for DNA, is basically achieved when the bsc0 α/γ refinement70 is applied. In the full set of Group III simulations (Table 3) considering different GQ folds, we observed the same dependence of BHB on ion type and parameters as was described for parallel GQs in the preceding part of this study (Supporting Information Figure S7). Thus, we suggest that all our above-discussed observations regarding BHB are valid for GQs with various folds. The effect of the GQ fold on the cation movements could not be confidently studied at the present simulation timescale, in line with the expected timescale of such processes.11 Ion exchanges in systems adherent to Group III were rare to provide reliable statistics. More specifically, we did not observe a single ion exchange in the Group III simulations except for simulations of 1KF1. It is not possible to tell whether it is a consequence of different GQ fold or of loops sterically hindering ion exchange. The stems of all the systems within Group III simulations were stable, but the tetrads experienced the usual BHB arrangements. In other words, except for some loop dynamics, the GQ stems in these simulations and on the present timescale 5−10 μs were entirely stable, albeit experiencing BHB. How To Select Cation Parameters for MD Simulations of G-Quadruplexes. Although we observed considerable dependence of analyzed GQ characteristics on the choice of ion type and parameters, it is not straightforward to recommend one parameter set that would be optimal for all GQ simulations. The choice of ion parameters may be affected by the investigated system and aims of a study. The model d[GGG]4 GQ was stable only with JC ions and partially with HFE KCl (let us remind that SPC/E water model specific ions are discussed). The other tested ions led to various disruptions of the system which we often considered as excessive. This supports the good performance of JC parameters in GQ simulations. While JC K+ ions did not exchange with the bulk, Na+ exhibited rapid exchange between the channel and bulk. Therefore, JC K+ parameters even enable one to study a system with different ions in the channel and bulk without exchanging them. Overall, we consider exchanges of JC Na+ perhaps too fast, and we would rather recommend the use of JC KCl for simulations of truncated systems, in combination with the SPC/E water model. Full GQs comprising loops were stable with all tested ions, and the rate of exchange of ions between the channel and the bulk was negligible. Thus, folded GQs can be simulated with diverse cation parameters. Still, intuitively, we would suggest that the use of JC parameters would be the safest option especially for long GQ simulations and should minimize the occurrence of some rare-event disruptions, since probably all force fields underestimate the stabilizing role of the cations on GQs due to the general approximation of the force field (see the Introduction).28 However, in studies aiming to investigate strand relocations or ion exchanges, one could consider utilization of other “less stable” cation parameters to speed up the dynamics. Nevertheless (despite that we lack rigorous experimental data about the rate of ion exchanges), we

which was stable during the simulation. IOD KCl and HFE KCl showed similar behavior − outer tetrads were in X-ray conformation and symmetric BHB prevailed for the inner tetrad. However, a few short flips back to the X-ray conformation were observed occasionally during the IOD KCl simulations. Simulations with DK KCl showed fluctuations between X-ray and BHB geometry of all tetrads. These flips were correlated, and the three tetrads were never in the X-ray conformation at the same time; if the inner tetrad flipped to the X-ray geometry, outer tetrads flipped to the bifurcated geometry and vice versa. Moreover, the 3′ tetrad also populated a conformation with bifurcation in the opposite direction, where N1 and N2 atoms of one guanine formed a bifurcated H-bond with O6 of a neighboring guanine (see Supporting Information Figure S8). Such a bifurcation leads to elongation of the distance between two diagonal O6 oxygens. Simulations of the truncated systems in NaCl showed the following trends: the 5′ end tetrad populated the symmetrically bifurcated geometry, the inner tetrad populated the asymmetrically bifurcated geometry, and the 3′ end tetrad populated either the X-ray geometry (in DNA simulations) or the symmetrically bifurcated geometry (in RNA simulations). Behavior of the inner tetrad in NaCl simulations was system specific. In full systems containing loops, asymmetric bifurcated geometry was stably locked in one of the two possible conformations. It should be noted that this behavior may be related to the presence of loops, which perturb the symmetry of the stems. However, frequent alterations of these conformations were observed in truncated systems. Asymmetry of bifurcation in truncated systems was more pronounced for IOD, HFE, and JC NaCl in comparison to AA NaCl, which is visible in Figure 3 as larger fluctuations of opposite O6 oxygen distances. Simulations of GQs with Different Folds Confirm the Same Ion Dependent Behavior of BHB as Was Observed for Parallel GQs. To verify the results, we have carried out the last set of computations (Group III, see Table 3), including several complete GQs with various folds, each simulated in three different force-field combinations: bsc0χOL4 DNA force field + JC NaCl, bsc0χOL4εζOL1 + JC NaCl, and OL15 + JC KCl. Note that all three DNA force-field versions share the basic AMBER Cornell et al. nucleic acid force field71 for the nonbonded terms and differ only in the degree of tuning of the backbone dihedral potentials. The latest OL15 includes reparametrization of all dihedral potentials compared to the Cornell et al. force field.50 Since the tetrad-ion interaction potentials were the same with all the DNA force-field versions and unchanged with respect to the original parametrization, it was meaningful to assume that the different DNA force-field versions should result in a similar description of cation interactions inside the stems. Full testing of this supposition was beyond our computational facilities. Nevertheless, for the full 1KF1 system with JC NaCl ion parameters we obtained trajectories with all three DNA forcefield versions. We did not observe any systematic effect of the DNA force field on analyzed parameters suggesting that the three backbone parametrizations affected neither the rate of ion exchanges nor the extent of presence of BHB. We nevertheless admit that the 1KF1 system contains only all anti-tetrads and canonical (B-DNA-like) backbone conformation. We cannot rule out that for some of the other topologies the details of the backbone parametrization might have some indirect effect on 3921


Article

Journal of Chemical Theory and Computation intuitively caution that “sampling enhancement” through intentional use of inaccurate force field can also lead to spurious structural changes. It has been noticed several times in the literature that simulations with the CHARMM force field show major instabilities of ions in the GQ stems,32 and this is the reason why we did not include this force field in our present tests. For a recent study revealing such ion expulsions see Ray et al.91 All the tested K+ parameters showed somewhat better geometries of tetrads in comparison to Na+, which led to asymmetric geometry with BHB. However, BHB was present to a certain extent in all simulations. Hence, if there is no clear preference for the type of salt, KCl provides better geometries of tetrads and less ion dynamics and could be a preferred option. Note that in the present study we did not test AA K+ parameters, since they have been shown previously to experience excessive ion expulsions in 0.5 μs simulations with the TIP3P water model and the bsc0 force field for the full 1KF1 GQ.42 In summary, the JC K+ parameters (at least with the SPC/E water model) represent an apparently safe option in studies of GQs. It should be noted that the JC parameters have a higher affinity to the DNA backbone compared to some other cation parametrizations.42,92 Such more subtle cation binding properties may have also some secondary effect on GQ simulations, for example on the dynamics of the loop regions.42 This is, however, beyond the scope of this study. Analysis of such more subtle (compared to the cation − G-tetrad coordination) effects may be complicated by the lack of unambiguous reference experimental data and difficulties of the force fields to describe the GQ loops.28,32,42,93−95 Effect of Water Model on the Presented Results. The details of GQ simulations may be also affected by the chosen water model to a certain extent.37,42 Several recent studies indeed revealed that the water model may influence MD simulations of nucleic acids.37,42,96−99 For example, it has been shown that the TIP5P model is not satisfactory for nucleic acids simulations with AMBER RNA force fields,97 while the OPC water model was shown to improve simulations of RNA tetranucleotides98 and calculations of stacking free energies in nicked DNA.99 Water models can affect stabilities and balance of various interactions in nucleic acids, such as base stacking and various types of hydrogen bonding. In addition, ions are usually parametrized to reproduce experimental properties like hydration free energies, ion−oxygen distance, and others.59 Thus, separate sets of parameters are commonly developed for different water models. As an example, JC parameters for K+ in SPC/E waters differ from JC parameters for K+ in TIP3P. Thus, ion parameters are often coupled with water models and both water model and ion parameters may affect GQ simulations. This complicates comparison of ion-dependent behavior of GQs in different water models. The present study primarily investigates the dynamics of desolvated ions inside GQ stems. Their influence is at first sight separated from the influence of the water molecules. However, the water model may affect other properties of the GQs, such as the groove widths, stacking stability, and movement of ions between the channel and bulk. To illustrate the potential sensitivity of the results to the water model and water-modelspecific ion parametrization, we performed a limited number of simulations in OPC77 and TIP3P78 water models with JC parameters for the d[GGG]4 system. We have first considered altogether three different parameter combinations with 10 × 5 μs of simulations in total (TIP4Pew JC ions58 were combined

with OPC). The full results are reported in the Supporting Information Figure S9 and the accompanying text. These simulations further confirmed the sensitivity of the d[GGG]4 system to the ion parameters. All three simulation sets indicated a lower stability of the d[GGG]4 stem than the fully stable T_1KF1_JC_KCl simulations. Roughly, in terms of structural stability, the KCl JC OPC, NaCl JC OPC, and KCl JC TIP3P simulations were similar to KCl IOD, NaCl AA, and KCl DK SPC/E simulations, respectively. The dependence of the BHB artifact and dynamics of ions on the ion type (Na+ vs K+) was consistent with the SPC/E simulations. Interesting dynamics was observed in two individual simulations (one in OPC and one in TIP3P water), where the 5′-tetrad refolded from an all anti to an all syn conformation. To the best of our knowledge, this is the first time that a complete rearrangement of the terminal tetrad has been spontaneously achieved in MD simulations, see the Supporting Information for details. To separate the effect of the ions inside the stem and of the water model, we further performed d[GGG]4 simulations with SPC/E specific JC parameters for K+ and Na+ in OPC water and vice versa (12 × 2.5 μs; see the Supporting Information). The simulations indicated a reduced structural stability of the d[GGG]4 in OPC water in comparison to SPC/E (see the Supporting Information for details). Although this effect was smaller compared to the effect of ion parametrizations, it was unambiguous. Dependence of BHB and ion exchange events on the type of ions (Na+ vs K+) was consistent with SPC/E simulations. The calculations thus suggest that water model and watermodel-specific parametrizations may further complicate the picture of structural dynamics of especially the d[GGG]4 system. Nevertheless, we suggest that most of the published literature data should not be dramatically affected by the water and ion models, since the simulations were done typically on completely folded GQs and on shorter timescales. Full systematic analysis of the role of the water model is beyond our current computer capabilities, as the number of possible models and their combinations with cation parameters is large. Nevertheless, any other combination of the parameters can be tested using the present d[GGG]4 data as a benchmark.

■

CONCLUDING REMARKS GQs are the most important noncanonical DNA molecules potentially occurring in telomeres and other parts of the genome. Likewise, RNA GQs can serve as regulatory elements in many biological processes. Because of their importance GQs are extensively studied by experimental and theoretical methods, including classical MD simulations. Desolvated monovalent cations present in the channel of GQs are necessary for stabilization of GQ stems through their interactions with G-tetrads. Monovalent ions are a key factor for folding and stability of GQs. It has been demonstrated in the past, that the overall electrostatic stabilization of GQs by the monovalent cations coordinated in the GQ stems is quite well described by contemporary MD simulations.28 However, finer details of the cation−GQ interactions are described less satisfactorily, causing some artifacts in the structural dynamics of the GQ structures, as in detail analyzed in our paper. Since there are various parametrizations of sodium and potassium ions available in the literature, we compared their performance for a set of GQ molecules, by analyzing 580 μs of simulation data, mostly in SPC/E water. Several different GQ 3922


Article


We could not test all the cation parametrizations available in the literature. As noted at the end of the Results and Discussion section, the comparison can be further complicated by considering different water models and coupling between ion parameters and water models. However, benchmarking of any other force field, ion and water parameters could be straightforwardly performed using the d[GGG]4 system and eventually extended to the r[ggg]4, 1KF1, and 3IBK systems. The analyses provided in the present paper can be used as the reference.

folds have been studied, and we assume that our results are generally valid for simulations of GQs. We have analyzed mainly the GQ structural stability, the rate of exchange of ions between the GQ channel and bulk solvent, and the presence of the bifurcated hydrogen bonding (BHB) artifact. All of the studied properties were highly dependent on the choice of ions (Figure 3). None of the studied ionic parametrizations was able to provide flawless performance in all analyzed GQ properties for all studied systems. We have made a number of specific observations that are detailed above. The identity and parametrization of ions strongly affected the behavior of the model d[GGG]4 GQ, which was fully stable only with JC SPC/E KCl ions and partially with HFE KCl (Figure 3 and Supporting Information Figure S3). In d[GGG]4 simulations, potassium ions did not show almost any exchanges between GQ channel and bulk, while sodium ions exchanged frequently (Figure 3). The exchanges depended also on the structural fluctuations that helped to cross the energy barrier for the cation exchange. Na+ ions tended to move inside the GQs in a concerted manner, while larger movements of the K+ ions were typically separated (Figure 8). The remaining studied RNA and DNA GQs were structurally stable (though suffering from the BHB artifact) regardless of the ion parametrization (Figure 5, Figure 6, Supporting Information Figures S4−S7). Even the parallel RNA GQ r[ggg]4 was structurally stable and stiffer than d[GGG]4. Thus, basic simulations of fully folded GQs can be safely done with a number of cation parametrizations, at least on the 10 μs timescale. Note that on the present simulation timescales, fully folded GQs are not expected to spontaneously lose ions, while exchanges of ions with the bulk should be rare.28 The potassium ion parametrizations showed somewhat better performance than sodium ions when considering the BHB artifact. Based on all data we suggest that the JC K+ parameters (at least with the SPC/E water model) represent an apparently safe option in studies of GQs. The reported sensitivity of the simulation data to the details of the cation parameters is quite surprising. Fortunately, as noted above, it does not significantly affect the most common simulations of fully folded cation-stabilized GQs. This is due to their structural stiffness and the large energy barrier that separates fully folded GQs from other parts of the folding landscape. Still, extensive sampling of the BHB G-tetrad arrangements is a reminder that the details of cation−tetrad interactions are captured imperfectly. This may have substantial consequences for other types of simulation studies, such as those monitoring details of the ion-binding processes and free energy calculations of the ion binding.47,100−103 Further, differences among various cation parameters may affect studies devoted to various aspects of GQ folding, when structures (and eventually transitions between them) other than the fully folded GQs are typically studied across the folding landscape.28 The choice of ions for a particular MD study could then be systemdependent and influenced by the aims of a given study. Preliminary investigation of the effect of the water model suggests that rate of ion exchanges and BHB artifact are not considerably affected by the water model, while the overall structural stability of the d[GGG]4 system is somewhat worsened in the OPC water model compared to the SPC/E. Water models can, for example, affect structural dynamics of the terminal G-tetrads exposed to the solvent. Thus, further systematic investigation of the effect of the water model on GQ simulations could be useful.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00257. Details of equilibration protocol; details of the MD simulation protocol; Group III GQs characterization; topologies of the simulated systems; description of development of selected d[GGG]4 simulations; graphs with development of selected properties of GQ in MD simulations of all systems; details on simulations with TIP3P and OPC water models (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Marek Havrila: 0000-0001-5066-8322 Barira Islam: 0000-0001-5882-6903 Michal Otyepka: 0000-0002-1066-5677 Funding

This work was supported by grant 16-13721S (M.H., B.I., P.S., and J.S.) and by the project SYMBIT reg. number: CZ.02.1.01/ 0.0/0.0/15_003/0000477 financed by the ERDF (J.S.). Further funding was obtained under the projects CEITEC 2020 (LQ1601, M.H. and J.S.) and LO1305 (P.S. and M.O.) of the Ministry of Education, Youth and Sports of the Czech Republic under the National Sustainability Programme II. J.S. acknowledges support by Praemium Academiae. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Sen, D.; Gilbert, W. Formation of Parallel Four-Stranded Complexes by Guanine-Rich Motifs in DNA and Its Implications for Meiosis. Nature 1988, 334, 364−366. (2) Sundquist, W. I.; Klug, A. Telomeric DNA Dimerizes by Formation of Guanine Tetrads between Hairpin Loops. Nature 1989, 342, 825−829. (3) Zahler, A. M.; Williamson, J. R.; Cech, T. R.; Prescott, D. M. Inhibition of Telomerase by G-Quartet DNA Structures. Nature 1991, 350, 718−720. (4) Siddiqui-Jain, A.; Grand, C. L.; Bearss, D. J.; Hurley, L. H. Direct Evidence for a G-Quadruplex in a Promoter Region and Its Targeting with a Small Molecule to Repress C-Myc Transcription. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 11593−11598. (5) Balasubramanian, S.; Hurley, L. H.; Neidle, S. Targeting GQuadruplexes in Gene Promoters: A Novel Anticancer Strategy? Nat. Rev. Drug Discovery 2011, 10, 261−275. (6) Rhodes, D.; Lipps, H. J. G-Quadruplexes and Their Regulatory Roles in Biology. Nucleic Acids Res. 2015, 43, 8627−8637. (7) Neidle, S. The Structures of Quadruplex Nucleic Acids and Their Drug Complexes. Curr. Opin. Struct. Biol. 2009, 19, 239−250.

3923


Article


Stacking Geometries: A Combined QM and MD Study. J. Phys. Chem. B 2013, 117, 9851−9856. (27) Yurenko, Y. P.; Novotny, J.; Mitoraj, M. P.; Sklenar, V.; Michalak, A.; Marek, R. Nucleic Acid Quadruplexes Based on 8-Halo9-Deazaxanthines: Energetics and Noncovalent Interactions in Quadruplex Stems. J. Chem. Theory Comput. 2014, 10, 5353−5365. (28) Šponer, J.; Bussi, G.; Stadlbauer, P.; Kührová, P.; Banás,̌ P.; Islam, B.; Haider, S.; Neidle, S.; Otyepka, M. Folding of Guanine Quadruplex Molecules−Funnel-Like Mechanism or Kinetic Partitioning? An Overview from MD Simulation Studies. Biochim. Biophys. Acta, Gen. Subj. 2017, 1861, 1246−1263. (29) Sponer, J.; Cang, X. H.; Cheatham, T. E. Molecular Dynamics Simulations of G-DNA and Perspectives on the Simulation of Nucleic Acid Structures. Methods 2012, 57, 25−39. (30) Spackova, N.; Berger, I.; Sponer, J. Nanosecond Molecular Dynamics Simulations of Parallel and Antiparallel Guanine Quadruplex DNA Molecules. J. Am. Chem. Soc. 1999, 121, 5519−5534. (31) Lane, A. N.; Chaires, J. B.; Gray, R. D.; Trent, J. O. Stability and Kinetics of G-Quadruplex Structures. Nucleic Acids Res. 2008, 36, 5482−5515. (32) Fadrna, E.; Spackova, N.; Sarzynska, J.; Koca, J.; Orozco, M.; Cheatham, T. E.; Kulinski, T.; Sponer, J. Single Stranded Loops of Quadruplex DNA as Key Benchmark for Testing Nucleic Acids Force Fields. J. Chem. Theory Comput. 2009, 5, 2514−2530. (33) Cieplak, P.; Dupradeau, F. Y.; Duan, Y.; Wang, J. Polarization Effects in Molecular Mechanical Force Fields. J. Phys.: Condens. Matter 2009, 21, 333102. (34) Baker, C. M. Polarizable Force Fields for Molecular Dynamics Simulations of Biomolecules. Wiley Interdisciplinary Reviews: Computational Molecular Science 2015, 5, 241−254. (35) Savelyev, A.; MacKerell, A. D., Jr. Balancing the Interactions of Ions, Water, and DNA in the Drude Polarizable Force Field. J. Phys. Chem. B 2014, 118, 6742−6757. (36) Islam, B.; Stadlbauer, P.; Krepl, M.; Koca, J.; Neidle, S.; Haider, S.; Sponer, J. Extended Molecular Dynamics of a c-Kit Promoter Quadruplex. Nucleic Acids Res. 2015, 43, 8673−8693. (37) Krepl, M.; Zgarbova, M.; Stadlbauer, P.; Otyepka, M.; Banas, P.; Koca, J.; Cheatham, T. E.; Jurecka, P.; Sponer, J. Reference Simulations of Noncanonical Nucleic Acids with Different Chi Variants of the Amber Force Field: Quadruplex DNA, Quadruplex RNA, and Z-DNA. J. Chem. Theory Comput. 2012, 8, 2506−2520. (38) Stadlbauer, P.; Krepl, M.; Cheatham, T. E.; Koca, J.; Sponer, J. Structural Dynamics of Possible Late-Stage Intermediates in Folding of Quadruplex DNA Studied by Molecular Simulations. Nucleic Acids Res. 2013, 41, 7128−7143. (39) Stadlbauer, P.; Kuhrova, P.; Banas, P.; Koca, J.; Bussi, G.; Trantirek, L.; Otyepka, M.; Sponer, J. Hairpins Participating in Folding of Human Telomeric Sequence Quadruplexes Studied by Standard and T-REMD Simulations. Nucleic Acids Res. 2015, 43, 9626−9644. (40) Zhang, X.; Xu, C. X.; Di Felice, R.; Sponer, J.; Islam, B.; Stadlbauer, P.; Ding, Y.; Mao, L.; Mao, Z. W.; Qin, P. Z. Conformations of Human Telomeric G-Quadruplex Studied Using a Nucleotide-Independent Nitroxide Label. Biochemistry 2016, 55, 360− 372. (41) Stadlbauer, P.; Mazzanti, L.; Cragnolini, T.; Wales, D. J.; Derreumaux, P.; Pasquali, S.; Sponer, J. Coarse-Grained Simulations Complemented by Atomistic Molecular Dynamics Provide New Insights into Folding and Unfolding of Human Telomeric GQuadruplexes. J. Chem. Theory Comput. 2016, 12, 6077−6097. (42) Rebic, M.; Laaksonen, A.; Sponer, J.; Ulicny, J.; Mocci, F. Molecular Dynamics Simulation Study of Parallel Telomeric DNA Quadruplexes at Different Ionic Strengths: Evaluation of Water and Ion Models. J. Phys. Chem. B 2016, 120, 7380−7391. (43) Zeng, X.; Zhang, L.; Xiao, X.; Jiang, Y.; Guo, Y.; Yu, X.; Pu, X.; Li, M. Unfolding Mechanism of Thrombin-Binding Aptamer Revealed by Molecular Dynamics Simulation and Markov State Model. Sci. Rep. 2016, 6, 24065.

(8) Largy, E.; Mergny, J.-L.; Gabelica, V. Role of Alkali Metal Ions in G-Quadruplex Nucleic Acid Structure and Stability. In The Alkali Metal Ions: Their Role for Life; Sigel, A., Sigel, H., Sigel, R. K. O., Eds.; Springer International Publishing: Cham, 2016; pp 203−258, DOI: 10.1007/978-3-319-21756-7_7. (9) Dingley, A. J.; Peterson, R. D.; Grzesiek, S.; Feigon, J. Characterization of the Cation and Temperature Dependence of DNA Quadruplex Hydrogen Bond Properties Using High-Resolution NMR. J. Am. Chem. Soc. 2005, 127, 14466−14472. (10) Sket, P.; Crnugelj, M.; Kozminski, W.; Plavec, J. 15NH4+ Ion Movement inside D(G4T4G4)2 G-Quadruplex Is Accelerated in the Presence of Smaller Na+ Ions. Org. Biomol. Chem. 2004, 2, 1970−1973. (11) Sket, P.; Virgilio, A.; Esposito, V.; Galeone, A.; Plavec, J. Strand Directionality Affects Cation Binding and Movement within Tetramolecular G-Quadruplexes. Nucleic Acids Res. 2012, 40, 11047−11057. (12) Davis, J. T. G-Quartets 40 Years Later: From 5′-GMP to Molecular Biology and Supramolecular Chemistry. Angew. Chem., Int. Ed. 2004, 43, 668−698. (13) Hud, N. V.; Smith, F. W.; Anet, F. A.; Feigon, J. The Selectivity for K+ Versus Na+ in DNA Quadruplexes Is Dominated by Relative Free Energies of Hydration: A Thermodynamic Analysis by 1H NMR. Biochemistry 1996, 35, 15383−15390. (14) Plavec, J. Metal Ion Coordination in G-Quadruplexes. In Metal Complex−DNA Interactions; John Wiley & Sons, Ltd.: 2009; pp 55− 93, DOI: 10.1002/9781444312089.ch3. (15) Creze, C.; Rinaldi, B.; Haser, R.; Bouvet, P.; Gouet, P. Structure of a D(TGGGGT) Quadruplex Crystallized in the Presence of Li+ Ions. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2007, 63, 682−688. (16) Phillips, K.; Dauter, Z.; Murchie, A. I.; Lilley, D. M.; Luisi, B. The Crystal Structure of a Parallel-Stranded Guanine Tetraplex at 0.95 Å Resolution. J. Mol. Biol. 1997, 273, 171−182. (17) Trajkovski, M.; Sket, P.; Plavec, J. Cation Localization and Movement within DNA Thrombin Binding Aptamer in Solution. Org. Biomol. Chem. 2009, 7, 4677−4684. (18) Meyer, M.; Steinke, T.; Brandl, M.; Sühnel, J. Density Functional Study of Guanine and Uracil Quartets and of Guanine Quartet/Metal Ion Complexes. J. Comput. Chem. 2001, 22, 109−124. (19) Song, J.; Ji, C.; Zhang, J. Z. H. The Critical Effect of Polarization on the Dynamical Structure of Guanine Quadruplex DNA. Phys. Chem. Chem. Phys. 2013, 15, 3846−3854. (20) Sponer, J.; Mladek, A.; Spackova, N.; Cang, X. H.; Cheatham, T. E.; Grimme, S. Relative Stability of Different DNA Guanine Quadruplex Stem Topologies Derived Using Large-Scale QuantumChemical Computations. J. Am. Chem. Soc. 2013, 135, 9785−9796. (21) Gkionis, K.; Kruse, H.; Platts, J. A.; Mladek, A.; Koca, J.; Sponer, J. Ion Binding to Quadruplex DNA Stems. Comparison of MM and QM Descriptions Reveals Sizable Polarization Effects Not Included in Contemporary Simulations. J. Chem. Theory Comput. 2014, 10, 1326− 1340. (22) Meyer, M.; Hocquet, A.; Suhnel, J. Interaction of Sodium and Potassium Ions with Sandwiched Cytosine-, Guanine-, Thymine-, and Uracil-Base Tetrads. J. Comput. Chem. 2005, 26, 352−364. (23) Woiczikowski, P. B.; Kubar, T.; Gutierrez, R.; Cuniberti, G.; Elstner, M. Structural Stability Versus Conformational Sampling in Biomolecular Systems: Why Is the Charge Transfer Efficiency in G4DNA Better Than in Double-Stranded DNA? J. Chem. Phys. 2010, 133, 035103. (24) Ho, J.; Newcomer, M. B.; Ragain, C. M.; Gascon, J. A.; Batista, E. R.; Loria, J. P.; Batista, V. S. Mod-QM/MM Structural Refinement Method: Characterization of Hydrogen Bonding in the Oxytricha Nova G-Quadruplex. J. Chem. Theory Comput. 2014, 10, 5125−5135. (25) Improta, R. Quantum Mechanical Calculations Unveil the Structure and Properties of the Absorbing and Emitting Excited Electronic States of Guanine Quadruplex. Chem. - Eur. J. 2014, 20, 8106−8115. (26) Lech, C. J.; Phan, A. T.; Michel-Beyerle, M. E.; Voityuk, A. A. Electron-Hole Transfer in G-Quadruplexes with Different Tetrad 3924


Article


(65) Lim, K. W.; Amrane, S.; Bouaziz, S.; Xu, W.; Mu, Y.; Patel, D. J.; Luu, K. N.; Phan, A. T. Structure of the Human Telomere in K+ Solution: A Stable Basket-Type G-Quadruplex with Only Two GTetrad Layers. J. Am. Chem. Soc. 2009, 131, 4301−4309. (66) Lim, K. W.; Alberti, P.; Guedin, A.; Lacroix, L.; Riou, J. F.; Royle, N. J.; Mergny, J. L.; Phan, A. T. Sequence Variant (CTAGGG)N in the Human Telomere Favors a G-Quadruplex Structure Containing a G.C.G.C Tetrad. Nucleic Acids Res. 2009, 37, 6239−6248. (67) Dai, J.; Punchihewa, C.; Ambrus, A.; Chen, D.; Jones, R. A.; Yang, D. Structure of the Intramolecular Human Telomeric GQuadruplex in Potassium Solution: A Novel Adenine Triple Formation. Nucleic Acids Res. 2007, 35, 2440−2450. (68) Dai, J.; Carver, M.; Punchihewa, C.; Jones, R. A.; Yang, D. Structure of the Hybrid-2 Type Intramolecular Human Telomeric GQuadruplex in K+ Solution: Insights into Structure Polymorphism of the Human Telomeric Sequence. Nucleic Acids Res. 2007, 35, 4927− 4940. (69) Lim, K. W.; Ng, V. C.; Martin-Pintado, N.; Heddi, B.; Phan, A. T. Structure of the Human Telomere in Na+ Solution: An Antiparallel (2 + 2) G-Quadruplex Scaffold Reveals Additional Diversity. Nucleic Acids Res. 2013, 41, 10556−10562. (70) Wei, D.; Parkinson, G. N.; Reszka, A. P.; Neidle, S. Crystal Structure of a c-Kit Promoter Quadruplex Reveals the Structural Role of Metal Ions and Water Molecules in Maintaining Loop Conformation. Nucleic Acids Res. 2012, 40, 4691−4700. (71) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A 2nd Generation Force-Field for the Simulation of Proteins, Nucleic-Acids, and Organic-Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (72) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham, T. E., III; Laughton, C. A.; Orozco, M. Refinement of the Amber Force Field for Nucleic Acids: Improving the Description of A/Γ Conformers. Biophys. J. 2007, 92, 3817−3829. (73) Zgarbová, M.; Otyepka, M.; Šponer, J. i.; Mládek, A. t.; Banás,̌ P.; Cheatham, T. E.; Jurečka, P. Refinement of the Cornell Et Al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. J. Chem. Theory Comput. 2011, 7, 2886−2902. (74) Zgarbová, M.; Luque, F. J.; Šponer, J.; Cheatham, T. E.; Otyepka, M.; Jurečka, P. Toward Improved Description of DNA Backbone: Revisiting Epsilon and Zeta Torsion Force Field Parameters. J. Chem. Theory Comput. 2013, 9, 2339−2354. (75) Galindo-Murillo, R.; Robertson, J. C.; Zgarbová, M.; Šponer, J.; Otyepka, M.; Jurečka, P.; Cheatham, T. E. Assessing the Current State of Amber Force Field Modifications for DNA. J. Chem. Theory Comput. 2016, 12, 4114−4127. (76) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (77) Izadi, S.; Anandakrishnan, R.; Onufriev, A. V. Building Water Models: A Different Approach. J. Phys. Chem. Lett. 2014, 5, 3863− 3871. (78) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (79) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (80) Hopkins, C. W.; Le Grand, S.; Walker, R. C.; Roitberg, A. E. Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning. J. Chem. Theory Comput. 2015, 11, 1864−1874. (81) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N· Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (82) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593.

(44) Li, H.; Cao, E. H.; Gisler, T. Force-Induced Unfolding of Human Telomeric G-Quadruplex: A Steered Molecular Dynamics Simulation Study. Biochem. Biophys. Res. Commun. 2009, 379, 70−75. (45) Buscaglia, R.; Miller, M. C.; Dean, W. L.; Gray, R. D.; Lane, A. N.; Trent, J. O.; Chaires, J. B. Polyethylene Glycol Binding Alters Human Telomere G-Quadruplex Structure by Conformational Selection. Nucleic Acids Res. 2013, 41, 7934−7946. (46) Le, H. T.; Dean, W. L.; Buscaglia, R.; Chaires, J. B.; Trent, J. O. An Investigation of G-Quadruplex Structural Polymorphism in the Human Telomere Using a Combined Approach of Hydrodynamic Bead Modeling and Molecular Dynamics Simulation. J. Phys. Chem. B 2014, 118, 5390−5405. (47) Bergues-Pupo, A. E.; Arias-Gonzalez, J. R.; Moron, M. C.; Fiasconaro, A.; Falo, F. Role of the Central Cations in the Mechanical Unfolding of DNA and RNA G-Quadruplexes. Nucleic Acids Res. 2015, 43, 7638−7647. (48) Stefl, R.; Cheatham, T. E.; Spackova, N.; Fadrna, E.; Berger, I.; Koca, J.; Sponer, J. Formation Pathways of a Guanine-Quadruplex DNA Revealed by Molecular Dynamics and Thermodynamic Analysis of the Substates. Biophys. J. 2003, 85, 1787−1804. (49) Rebic, M.; Mocci, F.; Laaksonen, A.; Ulicny, J. Multiscale Simulations of Human Telomeric G-Quadruplex DNA. J. Phys. Chem. B 2015, 119, 105−113. (50) Zgarbová, M.; Šponer, J.; Otyepka, M.; Cheatham, T. E.; Galindo-Murillo, R.; Jurečka, P. Refinement of the Sugar−Phosphate Backbone Torsion Beta for Amber Force Fields Improves the Description of Z- and B-DNA. J. Chem. Theory Comput. 2015, 11, 5723−5736. (51) Aznauryan, M.; Sondergaard, S.; Noer, S. L.; Schiott, B.; Birkedal, V. A Direct View of the Complex Multi-Pathway Folding of Telomeric G-Quadruplexes. Nucleic Acids Res. 2016, 44, 11024−11032. (52) Luo, D.; Mu, Y. Computational Insights into the Stability and Folding Pathways of Human Telomeric DNA G-Quadruplexes. J. Phys. Chem. B 2016, 120, 4912−4926. (53) Cavallari, M.; Garbesi, A.; Di Felice, R. Porphyrin Intercalation in G4-DNA Quadruplexes by Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113, 13152−13160. (54) Novotny, J.; Kulhanek, P.; Marek, R. Biocompatible XanthineQuadruplex Scaffold for Ion-Transporting DNA Channels. J. Phys. Chem. Lett. 2012, 3, 1788−1792. (55) Aaqvist, J. Ion-Water Interaction Potentials Derived from Free Energy Perturbation Simulations. J. Phys. Chem. 1990, 94, 8021−8024. (56) Smith, D. E.; Dang, L. X. Computer Simulations of NaCl Association in Polarizable Water. J. Chem. Phys. 1994, 100, 3757− 3766. (57) Dang, L. X.; Kollman, P. A. Free Energy of Association of the K+:18-Crown-6 Complex in Water: A New Molecular Dynamics Study. J. Phys. Chem. 1995, 99, 55−58. (58) Joung, I. S.; Cheatham, T. E. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020−9041. (59) Li, P.; Song, L. F.; Merz, K. M. Systematic Parameterization of Monovalent Ions Employing the Nonbonded Model. J. Chem. Theory Comput. 2015, 11, 1645−1657. (60) Case, D. A.; Babin, V.; Berryman, J. T.; Betz, R. M.; Cai, Q.; Cerutti, D. S.; Cheatham, T. E.; Darden, T. A.; Duke, R. E.; Gohlke, H., et al. Amber 14; University of California: San Francisco, CA, 2014. (61) Parkinson, G. N.; Lee, M. P. H.; Neidle, S. Crystal Structure of Parallel Quadruplexes from Human Telomeric DNA. Nature 2002, 417, 876−880. (62) Collie, G. W.; Haider, S. M.; Neidle, S.; Parkinson, G. N. A Crystallographic and Modelling Study of a Human Telomeric RNA (TERRA) Quadruplex. Nucleic Acids Res. 2010, 38, 5569−5580. (63) Parkinson, G. N.; Lee, M. P.; Neidle, S. Crystal Structure of Parallel Quadruplexes from Human Telomeric DNA. Nature 2002, 417, 876−880. (64) Wang, Y.; Patel, D. J. Solution Structure of the Human Telomeric Repeat D[AG3(T2AG3)3] G-Tetraplex. Structure 1993, 1, 263−282. 3925


Article

Journal of Chemical Theory and Computation (83) Havrila, M.; Zgarbová, M.; Jurečka, P.; Banás,̌ P.; Krepl, M.; Otyepka, M.; Šponer, J. Microsecond-Scale MD Simulations of HIV-1 DIS Kissing-Loop Complexes Predict Bulged-in Conformation of the Bulged Bases and Reveal Interesting Differences between Available Variants of the Amber Rna Force Fields. J. Phys. Chem. B 2015, 119, 15176−15190. (84) Cang, X. H.; Sponer, J.; Cheatham, T. E. Explaining the Varied Glycosidic Conformational, G-Tract Length and Sequence Preferences for Anti-Parallel G-Quadruplexes. Nucleic Acids Res. 2011, 39, 4499− 4512. (85) Islam, B.; Stadlbauer, P.; Neidle, S.; Haider, S.; Sponer, J. Can We Execute Reliable MM-PBSA Free Energy Computations of Relative Stabilities of Different Guanine Quadruplex Folds? J. Phys. Chem. B 2016, 120, 2899−2912. (86) Tong, X.; Lan, W.; Zhang, X.; Wu, H.; Liu, M.; Cao, C. Solution Structure of All Parallel G-Quadruplex Formed by the Oncogene RET Promoter Sequence. Nucleic Acids Res. 2011, 39, 6753−6763. (87) Clark, G. R.; Pytel, P. D.; Squire, C. J. The High-Resolution Crystal Structure of a Parallel Intermolecular DNA G-4 Quadruplex/ Drug Complex Employing Syn Glycosyl Linkages. Nucleic Acids Res. 2012, 40, 5731−5738. (88) Pagano, B.; Mattia, C. A.; Cavallo, L.; Uesugi, S.; Giancola, C.; Fraternali, F. Stability and Cations Coordination of DNA and RNA 14-Mer G-Quadruplexes: A Multiscale Computational Approach. J. Phys. Chem. B 2008, 112, 12115−12123. (89) Richardson, J. S.; Schneider, B.; Murray, L. W.; Kapral, G. J.; Immormino, R. M.; Headd, J. J.; Richardson, D. C.; Ham, D.; Hershkovits, E.; Williams, L. D.; et al. RNA Backbone: Consensus AllAngle Conformers and Modular String Nomenclature (an RNA Ontology Consortium Contribution). RNA 2008, 14, 465−481. (90) Bazzi, S.; Novotný, J.; Yurenko, Y. P.; Marek, R. Designing a New Class of Bases for Nucleic Acid Quadruplexes and QuadruplexActive Ligands. Chem. - Eur. J. 2015, 21, 9414−9425. (91) Ray, A.; Panigrahi, S.; Bhattacharyya, D. A Comparison of Four Different Conformations Adopted by Human Telomeric G-Quadruplex Using Computer Simulations. Biopolymers 2016, 105, 83−99. (92) Noy, A.; Soteras, I.; Luque, F. J.; Orozco, M. The Impact of Monovalent Ion Force Field Model in Nucleic Acids Simulations. Phys. Chem. Chem. Phys. 2009, 11, 10596−10607. (93) Islam, B.; Stadlbauer, P.; Krepl, M.; Koca, J.; Neidle, S.; Haider, S.; Sponer, J. Extended Molecular Dynamics of a c-Kit Promoter Quadruplex. Nucleic Acids Res. 2015, 43, 8673−8693. (94) Islam, B.; Sgobba, M.; Laughton, C.; Orozco, M.; Sponer, J.; Neidle, S.; Haider, S. Conformational Dynamics of the Human Propeller Telomeric DNA Quadruplex on a Microsecond Time Scale. Nucleic Acids Res. 2013, 41, 2723−2735. (95) Islam, B.; Stadlbauer, P.; Gil-Ley, A.; Pérez-Hernández, G.; Haider, S.; Neidle, S.; Bussi, G.; Banas, P.; Otyepka, M.; Sponer, J. Exploring the Dynamics of Propeller Loops in Human Telomeric DNA Quadruplexes Using Atomistic Simulations. J. Chem. Theory Comput. 2017, 13, 2458−2480. (96) Sklenovský, P.; Florová, P.; Banás,̌ P.; Réblová, K.; Lankaš, F.; Otyepka, M.; Šponer, J. Understanding RNA Flexibility Using Explicit Solvent Simulations: The Ribosomal and Group I Intron Reverse Kink-Turn Motifs. J. Chem. Theory Comput. 2011, 7, 2963−2980. (97) Kührová, P.; Otyepka, M.; Šponer, J.; Banás,̌ P. Are Waters around RNA More Than Just a Solvent? − an Insight from Molecular Dynamics Simulations. J. Chem. Theory Comput. 2014, 10, 401−411. (98) Bergonzo, C.; Cheatham, T. E. Improved Force Field Parameters Lead to a Better Description of RNA Structure. J. Chem. Theory Comput. 2015, 11, 3969−3972. (99) Hase, F.; Zacharias, M. Free Energy Analysis and Mechanism of Base Pair Stacking in Nicked DNA. Nucleic Acids Res. 2016, 44, 7100− 7108. (100) Akhshi, P.; Mosey, N. J.; Wu, G. Free-Energy Landscapes of Ion Movement through a G-Quadruplex DNA Channel. Angew. Chem., Int. Ed. 2012, 51, 2850−2854.

(101) Cavallari, M.; Calzolari, A.; Garbesi, A.; Di Felice, R. Stability and Migration of Metal Ions in G4-Wires by Molecular Dynamics Simulations. J. Phys. Chem. B 2006, 110, 26337−26348. (102) Akhshi, P.; Acton, G.; Wu, G. Molecular Dynamics Simulations to Provide New Insights into the Asymmetrical Ammonium Ion Movement inside of the [D(G3T4G4)]2 G-Quadruplex DNA Structure. J. Phys. Chem. B 2012, 116, 9363−9370. (103) Reshetnikov, R. V.; Sponer, J.; Rassokhina, O. I.; Kopylov, A. M.; Tsvetkov, P. O.; Makarov, A. A.; Golovin, A. V. Cation Binding to 15-TBA Quadruplex DNA Is a Multiple-Pathway Cation-Dependent Process. Nucleic Acids Res. 2011, 39, 9789−9802.

3926


Effect of Monovalent Ion Parameters on Molecular Dynamics

Recommend Documents