Derivation of Reliable Geometries in QM Calculations of DNA

Article pubs.acs.org/JCTC

Derivation of Reliable Geometries in QM Calculations of DNA Structures: Explicit Solvent QM/MM and Restrained Implicit Solvent QM Optimizations of G‑Quadruplexes Konstantinos Gkionis,† Holger Kruse,† and Jiří Šponer*,†,‡ †

Institute of Biophysics, Academy of Sciences of the Czech Republic, Královopolská 135, 612 65 Brno, Czech Republic CEITEC − Central European Institute of Technology, Masaryk University, Campus Bohunice, Kamenice 5, 625 00 Brno, Czech Republic

‡

S Supporting Information *

ABSTRACT: Modern dispersion-corrected DFT methods have made it possible to perform reliable QM studies on complete nucleic acid (NA) building blocks having hundreds of atoms. Such calculations, although still limited to investigations of potential energy surfaces, enhance the portfolio of computational methods applicable to NAs and offer considerably more accurate intrinsic descriptions of NAs than standard MM. However, in practice such calculations are hampered by the use of implicit solvent environments and truncation of the systems. Conventional QM optimizations are spoiled by spurious intramolecular interactions and severe structural deformations. Here we compare two approaches designed to suppress such artifacts: partially restrained continuum solvent QM and explicit solvent QM/MM optimizations. We report geometry relaxations of a set of diverse double-quartet guanine quadruplex (GQ) DNA stems. Both methods provide neat structures without major artifacts. However, each one also has distinct weaknesses. In restrained optimizations, all errors in the target geometries (i.e., low-resolution X-ray and NMR structures) are transferred to the optimized geometries. In QM/MM, the initial solvent configuration causes some heterogeneity in the geometries. Nevertheless, both approaches represent a decisive step forward compared to conventional optimizations. We refine earlier computations that revealed sizable differences in the relative energies of GQ stems computed with AMBER MM and QM. We also explore the dependence of the QM/MM results on the applied computational protocol.

■

INTRODUCTION Computational chemistry has been used extensively to study nucleic acid molecules, i.e. RNA and DNA. Most such studies conducted over the last 20 years have relied on atomistic explicit solvent simulations utilizing pair-additive molecular mechanics (MM) force fields.1,2 More expensive quantum chemical (QM) calculations have been used to study small model systems such as base pairs, base stacks, ion-binding patterns, and backbone conformations.3−38 QM calculations also played an essential role in tuning nucleic acid force fields.39−42 An important advantage of QM methods is that they make it possible to study effects that cannot be described by classical MM.11 Because of this and the increasingly apparent limitations of force field-based methods,43 there is a growing interest in applying QM methods to larger nucleic acid fragments. This has been facilitated by ongoing improvements in computational quantum chemistry methods and the performance of commodity computing hardware that have made it possible to apply such sophisticated approaches to ever larger systems. QM calculations cannot replace extended MM-based explicit solvent simulations of nucleic acids because they are unsuitable for sampling conformational space, which is absolutely essential in biomolecular modeling. However, they © 2016 American Chemical Society

can considerably improve the accuracy of predictions based on MM calculations and facilitate further tuning of force fields. Several recent publications44−46 have reported geometry optimizations and energy calculations based on the application of modern dispersion-corrected DFT to essentially complete nucleic acid building blocks including two-quartet guanine quadruplex (GQ) stems,44 three-base pair B-DNA helices,45 and the 11-nucleotide Sarcin-Ricin RNA motif.46 By complete nucleic acid building blocks, we mean systems that are large enough to autonomously adopt stable biochemically relevant three-dimensional configurations under real experimental conditions; in practice, this means systems of at least 6−10 nucleotides. Modern dispersion-corrected DFT methods (i.e., methods that include London dispersion interactions within the framework of Kohn−Sham DFT) achieve accuracies that undoubtedly surpass those of pair-additive MM, which is not the case for alternatives such as semiempirical QM calculations. Note that traditional DFT methods do not properly account for the dispersion energy, which has been a major problem in biochemically oriented applications of DFT methods for some Received: October 30, 2015 Published: February 25, 2016 2000

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation time. At present the best-established dispersion-corrected approaches are probably DFT-D347,48 and the Minnesota family of functionals.49,50 The former approach is used in this work. These and other DFT dispersion corrections have been reviewed by Grimme,51 and extensive benchmarking studies on dispersion-corrected DFT have been presented by Goerigk et al.52−54 and others.55,56 Data from DFT-D3 calculations on GQ stems44 have been used to correct earlier MD-based free energy calculations on the relative stability of GQ stems,57 leading to greatly improved agreement with experimental data. Recent studies have demonstrated that DFT calculations on complete nucleic acid building blocks are technically feasible. However, they also showed that it is not straightforward to achieve biochemical relevance with QM computations. Standard QM optimizations often yield geometries that deviate significantly from those that are biochemically relevant because the studied systems are typically incomplete, having implicit solvent environments and being without surrounding nucleic acid elements.44−46,58 This leads to artifacts such as non-native H-bonds and spurious backbone conformations (see below for more details). Such unrealistic geometries are entirely outside the conformational space sampled by real molecules and if they are used in subsequent calculations of energies and other properties, the transferability of the results will be severely limited. The difficulty of obtaining realistic geometries is currently the main factor limiting the routine use of large-scale QM computations in the study of nucleic acids. We recently introduced an efficient method that reduces uncertainties in QM calculations on biomolecular building blocks, namely restrained optimizations. 59 By applying harmonic restraining potentials to geometric parameters such as the backbone dihedrals of nucleic acids, we can gently relax the biomolecular building blocks while avoiding non-native structural features and artifacts. Conventional constraints are unsuitable for such optimizations because some residual flexibility is required when studying large systems. Here, we compare this method to an alternative - hybrid QM/MM computations - by performing unrestrained QM/MM optimizations with explicit solvation and restrained QM optimizations with implicit solvation for different topologies of selected two-quartet GQ DNA stems. GQ stems represent ideal systems for testing these methodologies due to their unique combination of topological variability and stiffness.43,44 G-quadruplexes are the most important class of noncanonical DNA molecules. Guanine-rich DNA regions are associated with the formation of units of planar guanine quartets that stack to form GQ stems. The stacked arrangement of the quartets is characterized by a central channel that hosts stem-stabilizing monovalent cations (K+, Na+, or NH4+). The stem can be formed by the intermolecular combination of different DNA strands or intramolecular folding of a single strand. In the latter case, the G-stretches forming the stems can be connected by diagonal, lateral, or propeller loops, giving rise to diverse architectures. The type of the loop dictates the stem’s overall topology and is inter-related with the characteristic permitted alternations of anti and syn conformations of the constituent guanosines. Together these factors give rise to a multitude of possible GQ topologies,60−67 and this topological variability is crucial in their biological functions. An example of a GQ stem is shown in Figure 1. GQs can form in telomeric regions of the genome and during biological processes such as transcription and replication.68−70

Figure 1. Example of an antiparallel two-quartet GQ stem with syn and anti orientations (defined by the rotation of the base around the sugarbase bond) of the guanines with respect to the sugar rings. Black arrows indicate the 5′-3′ directionality of the four strands.

Additionally, GQs can serve as aptamers71,72 and potential drug targets61,73,74 or even be used in molecular assemblies relevant to nanotechnology applications.75−77 Because of their biological importance and potential practical applications, the stability, variability, and folding mechanisms of GQ DNA have been studied intensely.60−63,78−85 Biologically relevant GQs are also readily formed by RNA molecules.86−88 Many experimental methods have been used to investigate GQs, such as NMR, X-ray crystallography, Circular Dichroism, Mass Spectrometry, and other methods.64,89,90 For some time, molecular dynamics (MD) simulations have been the main computational tool used to study aspects of GQs such as their folding, stability, and ion dynamics.43,57,91−98 However, it is becoming common for such simulations to be performed with time scales of microseconds, making the limitations of existing force fields increasingly apparent and justifying the use of QM computations wherever possible. QM calculations were used to study the properties of guanine quartets and ion binding to the quartets.8,11−22 As noted above, our recent DFT-D3 study on complete twoquartet GQ stems (including the sugar−phosphate backbone) revealed, in combination with MD simulation data, the relative free energies of several different GQ topologies.44 The present study substantially extends this work and demonstrates how GQ stems can be investigated in a highly controllable manner with the aid of restrained continuum solvent QM optimizations and explicit solvent QM/MM optimizations. The two methods have distinct advantages and limitations, but both provide QMrelaxed geometries that are free of non-native artifacts. The continuum-solvent restrained optimizations typically yield optimized structures that are close to the target geometries. However, this is also a weakness of the method because we often lack unambiguous target structures due to the genuine flexibility of nucleic acids and uncertainty in the experimental structures caused by resolution limits and refinement errors. QM/MM optimizations are performed in a realistic explicit solvent environment and, at least for GQ stems, do not require any restraints. However, the specific starting configuration chosen for the solvent can produce random heterogeneities in the optimized structures. This problem is not alleviated even when MD simulations are performed to equilibrate the water molecules because the QM/MM optimization is still initiated from a single explicit solvent and solute configuration. In addition, proper relaxation of the solvent would have to be coupled with some concerted relaxation of the DNA molecule, i.e. of the QM region. It is important to be aware that real 2001

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation

Figure 2. Schematic representations of the two-quartet GQ DNA topologies considered in this work. The known multiple-quartet quadruplex topologies can be constructed using these two-quartet stems as the basic building blocks. Guanine nucleosides are represented by rectangles (white for anti, gray for syn), and arrows show the 5′-3′ directionality of each strand. Note that two distinct starting configurations are investigated for the SA-abab fold.

in implicit solvent environments using the Poisson−Boltzmann surface area approach). This MM-based free energy order was two years later corrected based on the differences between the QM and MM potential energies calculated for a series of optimized structures of two-quartet stems in continuum solvent - for full details see ref 44. The QM-based correction was essential to obtain agreement with the experimentally known preferences of the GQ quadruplexes to form AA stems when not restricted by single-stranded loops. The present work is a major extension of the latter study.44 Solvated Structures of Stems. Explicit water molecules and counterions were added using the LEaP program105 of the AMBER14 package.106 Water molecules were added within an octahedral box with closeness and distance parameters of 1.0 and 3.5 Å, respectively (see Table 1 for the number of waters added in each case). One counterion (either K+ or Na+, chosen to match the channel ion) per phosphate group was also added. The force field parameters used in the calculations are specified below. Thus, the total charge of the complete models was +1 (four MM counterions, one QM channel ion, and four phosphate groups). The positions of the counterions were those assigned automatically by the program. We added four rather than three “bulk” ions to avoid introducing nonuniformity due to the distribution of the bulk ions. For some QM/MM optimizations on the AA system, short MD runs with the solute frozen were performed to better equilibrate the solvent molecules (for more details see the Supporting Information). QM/MM Calculations. All QM/MM geometry optimizations were performed using the ONIOM method.107−110 In the simplest form of the ONIOM method the system consists of two layers (two-layer ONIOM): the high layer is described by a higher level method and the low layer by a lower level method. Unless otherwise specified, the high layer (QM) region consisted of the GQ stem and the channel ion and was treated at the DFT level using the TPSS111 functional with the TZVP112 basis set and D3 dispersion correction with the BJ damping function47,48 (TPSS-D3(BJ)/TZVP). The density fitting approximation was used to speed up the calculations. The low layer consisted of the water molecules and the bulk counterions and was treated at the MM level using the TIP3P113 water model and compatible ion parameters114 unless otherwise stated. The two-layer ONIOM method requires three separate calculations: QM and MM calculations on the high layer and an MM calculation on the whole system. In our case, this meant performing TPSS-D3(BJ)/TZVP and MM calculations on the stem and an MM calculation on the entire system. The total ONIOM energy (although not the main result of our work) is given by the extrapolation

molecules always sample a broad range of configurations during the genuine thermal fluctuations of the whole system. Thus, any calculation based on evaluation of a small number of configurations will necessarily suffer from some sampling limitations.

■

METHODS Initial Structures of GQ Stems. The input structures of the eight stems including the K+ channel ion are taken from ref 44, where more details about the structures can be found (Figure 2). We maintain the same model names as in the original paper. Syn and anti guanines are abbreviated as “S” and “A”, respectively, and eight models were considered: one model with four anti-anti steps (abbreviated as model AA), one model with four anti-syn steps (AS), five models with four syn-anti steps (SA-aaaa, SA-aaab, SA-aabb, SA-abab, SA-abab-2), and one model with three anti-anti steps and one syn-syn step (3AA +1SS), cf. Figure 2 and Table 1. Sequences of four lowercase Table 1. Details of the Structures Considered in This Work PDBa

exp.b

AA

352D99

X-ray

AS SA-aaaa SA-aaab

1JPQ100 3TVB101 2GKU102

X-ray X-ray NMR

SA-aabb SA-abab SA-abab-2 3AA+1SS

1JPQ100 148D103 2AVH104 2GKU102

X-ray NMR X-ray NMR

model

model constructionc 2nd, 3rd quartets of the 1st molecule 2nd, 3rd quartets 1st, 2nd 5′-end quartets 1st, 2nd quartets of the 1st model 1st, 2nd quartets 1st model 1st, 2nd quartets 2nd, 3rd quartets of the 1st model

no. of watersd 604 665 558 725 667 710 777 624

a

PDB code of the parent experimental structure. bType of experimental structure. cDetails of model construction from the experimental structure (taken from ref 44, see the Supporting Information of ref 44 for the coordinates). dNumber of explicit water molecules included in our QM/MM calculations.

letters are used to indicate the relative orientations of the strands in the different model topologies: strands oriented “upward” in Figure 2 are denoted “a” and those oriented “downward” are denoted “b”. In each case, the first letter of the sequence indicates the orientation of the strand in the front-left corner of the schematic representations shown in Figure 2, and the sequence proceeds in anticlockwise order. Structures with Na+ as the channel ion were prepared by replacing K+ ions in the original input geometries with Na+.44 A previous publication57 described the derivation of relative free energies for these two-quartet stems using explicit solvent MD simulations followed by MM-PBSA computations57 (free energy analysis of the MD trajectory based on MM calculations 2002

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation MM MM QM E ONIOM = E(stem) + E(stem + water + counterions) − E(stem)

DNA in the AMBER program package. For details of these parameters, see the AmberTools15 May 2015 update or http:// fch.upol.cz/ff_ol/. It should be noted that the AMBER code sometimes internally refers to nucleic acid force fields as being part of the ffxxSB biomolecular force field set. The ffxxSB abbreviation can be very confusing and ambiguous when used in the literature for nucleic acid force fields because it is primarily used to denote a set of unrelated protein force fields. Implicit solvent environment was modeled using the Poisson−Boltzmann method (PB). The PB optimizations (MMPB) were performed using the xopt code59 coupled to the AMBER14 package.106 The grid spacing for the numerical PB calculations was set to 0.2 Å, the arc resolution to 0.0625 Å, and the solvent probe radius to 1.6 Å. Calculations of Energies. The relative energies of different stem topologies were evaluated at their QM/MM optimized geometries, stripped of the water molecules and counterions (the bound channel ion was always retained in the structure) and at the restrained optimized geometries (hereafter abbreviated as “QMrest” geometries). The energies were calculated at the TPSS-D3(BJ)/def2-TZVP level using the COSMO implicit solvation model and the outlying charge correction (a correction for charge balancing).124 Additional calculations were performed using the SMD125 solvation model and the TZVP and def2-TZVP basis sets. SMD and COSMO calculations were performed using the Gaussian09 (Revision D.01)118 and ORCA (version 2.9.1)126 packages, respectively. Justification of the Choice of the TPSS-D3 Method. The meta-GGA TPSS-D3 is one of the most widely tested medium-cost dispersion-corrected DFT methods. In studies on nucleic acid systems, TPSS-D3 achieved results comparable to the more demanding global hybrid PW6B95-D3 method while having the advantage of permitting the structural relaxation of such systems due to its lower computational cost.44,46 Some more general benchmarking studies on TPSS-D3 are cited in the Introduction. We also note that hybrid functionals such as B3LYP are prohibitively expensive for geometry optimizations of systems exceeding 200 atoms, unless a small double-ζ basis set is used (e.g., B3LYP/6-31G*). However, the use of small basis sets would introduce significant intra- and intermolecular basis set superposition error, affecting both interaction energies and structures.127 MD Simulations. As explained in the Results and Discussion, we performed one classical MD simulation on the AA topology. This simulation was carried out using the protocol of Cang et al.57 and is summarized in the Supporting Information, under the heading “Protocol for the MD simulation on the AA-K+ topology”. Analysis of Results. The results were analyzed using our own scripts, the “cpptraj”128 utility of the AMBER Tools package, and the VMD129 program. Images were generated using VMD and the UCSF Chimera package.130 List of Abbreviations. To help the reader follow the discussion of the results, a list of abbreviations used in this study is provided in the Supporting Information.

(1)

Note that the stem always includes its internally bound QMtreated channel ion. The full MM calculations were performed using the parmbsc040 DNA parameters of the Cornell et al. force field.115 The channel ions were assigned the same parameters as the bulk ions. Our QM partitioning does not involve breaking of covalent bonds between the two layers. Therefore, the only interactions between the two layers are nonbonded and were accounted for at the MM level using a mechanical embedding (ME) scheme. We are aware that the description of the electronic structure of the negatively charged nucleic acids by the ME scheme may be imperfect if all the negative charge is included in the DFT region (as was the case in this work). However, the results obtained using this method were sufficiently accurate for our purposes here (see below). A more realistic approach to the problem would be to embed the point charges of the MM atoms into the QM Hamiltonian. Then the electrostatic interactions between the low and high layers would be accounted for at the QM level by the polarized Hamiltonian (electrostatic embedding, EE). However, EE optimizations are usually much harder to converge (see below), which is the main reason we used the ME approach in most of our computations. ME and EE geometry optimizations were performed for a number of test GQ stem systems, and in no case did the two schemes yield results that differed significantly. We therefore assume that all the ME-optimized geometries for the systems considered in this work are acceptably accurate. Microiterations116 were used in all geometry optimizations, and the trust radius was not updated at each optimization step, in keeping with earlier work.117 All ONIOM QM/MM calculations were performed using the Gaussian09 (Revision D.01) software.118 Restrained QM Optimizations. We also performed restrained optimizations at the TPSS-D3(BJ)/TZVP119 level with the COSMO120−122 continuum solvation model for water (ε = 78.4).59 The backbone angles (α−ζ and χ) of all four strands were restrained to their input values using a harmonic penalty function with a restraint constant of 0.01 Eh/rad2. The optimizations were performed using the open source code xopt59 for energy and gradient calculations, which was interfaced with Turbomole 6.3 as an external optimizer.123 The Turbomole grid m4 was used for the quadrature integration. The convergence threshold for the SCF energy and the energy change of successive optimization cycles was 10−7 Eh; that for the gradient norm was 10−3. Force Field Optimizations. We performed MM optimizations using the 2007 parmbsc040 reparametrization of the basic Cornell et al.115 ff99 force field, upgraded with the 2012 χOL441 and 2013 εζOL139 refinements. The bsc0 parameters provide improved α/γ dihedral potentials to stabilize DNA simulations, the χOL4 parameters provide optimized χ dihedrals to improve simulations of GQs and other DNAs with syn nucleotides, and the εζOL1 parameters provide tuned ε/ζ dihedral potentials that correct the DNA helical twist and populations of BI/BII backbone substates. For a recent overview of these force field modifications, see ref 2. Since all of the dihedral refinements are one-dimensional (uncoupled) tweaks of the backbone dihedral potentials, they are separable and, in principle, can be used in different combinations. However, the bsc0χOL4εζOL1 combination featuring all three sets of refinements is currently the recommended version for

■

RESULTS AND DISCUSSION Earlier MM and QM Computations on Two-Quartet GQ Stems. Before discussing the results of our new calculations, we will briefly review previously reported calculations on the same GQ stems because their limitations motivated the present work.44,57 Each two-quartet GQ stem consists of four GpG steps (dinucleotides). Estimates of the 2003

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation

Figure 3. Backbone-base H-bonds observed in our computations. Left: syn conformation displaying an O5′-H···G(N3) interaction. Middle: anti conformation with the spurious N2−H···O4′ interaction. Right: anti conformation showing two weak C8−H···O5′ interactions.

non-native conformations during QM optimizations of DNA building blocks: unconstrained QM/MM explicit solvent optimizations and restrained implicit solvent QM optimizations. These methods and the corresponding geometries are henceforth denoted QM/MM and QMrest, respectively. We also explore the successive application of both methods, i.e. QM/ MM geometry optimization followed by restrained QM optimization and vice versa, hereafter denoted as QM/MM → QMrest and QMrest → QM/MM, respectively. Our results are consistent with the free energy predictions from our earlier work but have much lower levels of uncertainty than the previous calculations did. In addition, we provide some guidelines for performing QM optimizations of nucleic acid building blocks and verifying that the optimized structures are sufficiently accurate to perform additional analyses. Specific Intramolecular H-Bonds in GQ Stems. Before introducing the results, we will discuss the kinds of intramolecular H-bonds that may occur in the studied GQ structures and which affect the calculated energies of different topologies (Figure 3). The O5′-H···G(N3) interaction occurs exclusively in 5′-syn terminal guanosines. It has been estimated that one such Hbond stabilizes the structure by ∼7 kcal/mol in continuum solvent potential energy QM descriptions.44 The O5′-H···G(N3) H-bond exists in experimental structures when the guanine in question is the 5′-terminal residuum in the DNA strand possessing the free 5′-OH terminal group; X-ray structures reveal that the O5′···G(N3) distance in such cases is ∼3 Å. It is this interaction that forces the 5′-terminal guanine to adopt the syn conformation.44,57 Among the experimental structures used as initial structures in this work, short O5′···G(N3) distances are observed in all four strands of the SA-aaaa structure and one strand each of the SA-aabb and SA-abab-2 structures. In our QMrest optimizations, a short initial O5′···G(N3) distance invariably led to the formation of such H-bonds. In contrast, this H-bond does not form in the QM/ MM optimizations because we initially oriented the corresponding hydrogen away from the acceptor. However, the Hbond is observed in the structures arising from the QMrest → QM/MM optimizations (in which it persists throughout the optimization process), and it also occurs in a majority of the relevant nucleotides in the QM/MM → QMrest structures. The presence or absence of this interaction is noted explicitly when presenting the results because it affects the calculated relative energies (for further details, see ref 44). It is worth repeating that structures featuring this interaction and those that lack it are both experimentally relevant (see above). Thus, we did not attempt to unify their occurrence. A case in point is the SA-abab model derived from the NMR structure of the double-quartet thrombin binding aptamer

relative stabilities of 5′-anti-anti-3′ and 5′-syn-anti-3′ GpG steps embedded in different GQ stems obtained from earlier DFTD3 calculations44 differed very significantly from the corresponding relative stabilities obtained from equivalent MM calculations. The difference was attributed to the inaccuracy of the DNA force field and was most apparent when comparing QM and MM potential energies for the AA and four SA GQ stem structures (cf. Figure 2). These calculations explained a major discrepancy between an earlier MD-derived MM-PBSA free energy ordering of the different GQ stem topologies57 and the available experimental data. Specifically, the MM-PBSA calculations predicted the 5′-syn-anti-3′ GpG step to be more stable than the 5′-anti-anti-3′ GpG step by 3.5 kcal/mol.57 This MM-PBSA result was inconsistent with two well-established experimental findings, namely, the known propensity of GQs to form all-parallel all-anti stems and the fact that the 5′-anti-anti3′ and 5′-syn-anti-3′ GpG steps should be very close in energy. Thus, according to experiment, DNA GQs maximize the number of these two dinucleotide conformers (see the discussion of the experimental data in ref 44). Adding the difference between the DFT-D3 and MM potential energy calculations to the MM-PBSA free energy data as a corrective term clearly improved the agreement between the computations and experiments: with this correction applied, the 5′-antianti-3′ and 5′-syn-anti-3′ GpG steps were predicted to be isoenergetic within the calculations’ error margins.44 These results clearly demonstrated both the limitations on the accuracy of MM-based descriptions of nucleic acids and the utility of QM approaches. Despite this success, the previously reported continuum solvent QM-optimized geometries were not fully satisfactory, with severe deformations of some structures due to the formation of non-native (gas-phase like) intramolecular Hbonds and sugar−phosphate backbone deformations.44 This prompted us to perform additional QM optimizations with the inclusion of a few specific water molecules disrupting the nonnative interactions, as well as QM single-point energy calculations using artificially rigidified MM-optimized geometries.44 Although the energies obtained with these different sets of calculations were qualitatively consistent, the quantitative differences between them were sizable. The energy correction suggested in the preceding study (i.e., the difference between the QM and MM energies that was added to the MMPBSA free energy predictions for the relative stability of different GpG dinucleotides from ref 57) was thus based on ad hoc averaging of two very different sets of computations.44 The ad hoc nature of the corrections used to eliminate nonnative structural features in our preceding work was unsatisfactory. Therefore, in this work, we compare two systematic approaches intended to prevent the formation of 2004

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation Table 2. Relative Energies (ΔE) of Different G-Stem Topologies with Respect to the AA Topologyg K+ ΔEQM geometry

AA

AS

QMunresta QMunrestb QMunrest*c QM/MM QMrest QM/MM→QMrest QMrest→QM/MM

0.0 0.0 0.0 0.0 0.0 0.0 0.0

+16.7 +14.7 +8.4 +14.4 +12.0 +6.2 +17.9

bsc0

AA

SA-aaaa −21.5 −23.4 −29.0 +20.0 −24.9 −11.3 −26.5

[4] [4] [4] [4] [2] [4]

SA-aaab

SA-aabb

SA-abab

SA-abab-2

+3.9 +0.7 −4.9 +5.5 +0.8 −4.6 +3.0

+0.4 [1] −2.8 [1] −8.5 [1] +7.7 −2.4 [1] −7.7 [1] −3.9 [1] ΔEMMe

+9.8[1]d +7.4[1]d +1.8[1]d −3.0 [1] +5.1[1]d −7.1[1] +7.9[1]d

−4.4 −7.1 −12.7 +10.7 −8.5 −8.8 +2.9

[1] [1] [1] [1] [1] [1]

3AA+1SS +6.8 +5.6 0.0 +12.5 +14.5 +8.1 +14.7

AS

SA-aaaa

SA-aaab

SA-aabb

SA-abab

SA-abab-2

3AA+1SS

−1.5[1] −4.0 [1] −2.1[1]d −4.2 [1] −4.4[1]d SA-abab

−11.3 [1] −7.3 −13.9 [1] −14.5 [1] −13.4 [1] SA-abab-2

+11.1 +12.7 +13.0 +12.4 +12.8 3AA+1SS

−1.3[1]d −7.5 [1] −2.8 −6.9 [1] −5.7 SA-abab

−14.3 [1] −10.3 −14.3 [1] −16.3 [1] −14.6 [1] SA-abab-2

+18.2 +13.9 +16.5 +15.0 +16.1 3AA+1SS

QMunrest→MMPB QM/MM→MMPB QMrest→MMPB QM/MM→QMrest→MMPB QMrest→QM/MM→MMPB bsc0χOL4

0.0 0.0 0.0 0.0 0.0 AA

+19.6 +14.4 +18.0 +16.2 +19.5 AS

−26.2 [4] −2.5 −27.4 [4] −16.3 [2] −28.7 [4] SA-aaaa

−6.7 −7.6 −7.5 −7.6 −6.4 SA-aaab

−5.6 [1] −7.3 [1] −7.0 [1] −7.7 [1] −5.7 [1] SA-aabb

QMunrest→MMPB QM/MM→MMPB QMrest→MMPB QM/MM→QMrest→MMPB QMrest→QM/MM→MMPB bsc0χOL4εζOL1

0.0 0.0 0.0 0.0 0.0 AA

+15.3 +11.4 +16.1 +12.0 +16.5 AS

−25.6 [4] −4.5 −26.3 [4] −17.1 [2] −29.7 [4] SA-aaaa

−9.0 −10.8 −9.5 −10.5 −9.2 SA-aaab

−7.4 [1] −10.8 [1] −7.4 [1] −10.0 [1] −7.7 [1] SA-aabb

QMunrest→MMPB QM/MM→MMPB QMrest→MMopt QM/MM→QMrest→MMPB QMrest→QM/MM→MMPB earlier MMPBf

0.0 0.0 0.0 0.0 0.0 0.0

+15.1 +10.4 +13.8 +10.9 +15.0 −2.7

−24.5 −4.3 −27.0 −17.4 −29.4 −53.9

−7.7 −10.2 −9.8 −10.5 −8.8 −22.9

[4] [4] [2] [4] [4]

−7.0 −9.4 −7.7 −10.3 −7.7 −27.0

[1] [1] [1] [1] [1] [1]

d

0.0[1]d −7.1 [1] −4.6 −6.9 [1] −5.9 −14.8

−12.8 −10.2 −14.6 −16.4 −14.2 −31.3

[1] [1] [1] [1] [1]

+16.7 +12.9 +13.6 +13.1 +14.7 +2.1

Na+ ΔEQM QM/MM QMrest

AA

AS

SA-aaaa

SA-aaab

SA-aabb

SA-abab

SA-abab-2

3AA+1SS

0.0 0.0

+9.7 +12.1

+8.0 [1] −24.5 [4]

+4.1 +1.2

+7.6 −1.8 [1]

+3.5 +5.7[1]d

+8.8 −8.7 [1]

+15.0 +13.8

a TZVP data from ref 44. bdef2-TZVP (recalculated in this work). cΔE with respect to AA QMunrest* geometry. dAn unusual N2−H···O5′ H-bond is formed instead of the O5′-H···G(N3) interaction; see the text for an explanation. ePB ΔEs of MMPB-optimized geometries using the bsc0, bsc0χOL4, and bsc0χOL4εζOL1 force fields. The MMPB optimizations started from different QM geometries, as indicated. fA previously reported MMPBoptimized geometry derived with the bsc0 force field starting from the QMunrest geometry; these data were taken directly from ref 44. Note that the PB optimizations in ref 44 were done using the AMBER code with the version of the PB optimizer that was available at the time, which was very inefficient. The optimizations resulted essentially in mere relaxation of the bond lengths and bond angles, without any visible geometrical changes. For a full discussion see ref 44. gAll energies are calculated using the DFT-D3 method in COSMO solvent, mostly with the def2-TZVP basis set. QMunrest and QMunrest* are unconstrained geometries taken from ref 44, where * indicates optimizations with four explicit water molecules for the AA and AS structures. The MM data (MMPB) are obtained with the PB solvent and using the specified force field versions. Numbers in square brackets denote the number of O5′-H···G(N3) H-bonds present in the structure; note that one such H-bond confers approximately 7 kcal/mol of stabilization. Coordinates for all optimized structures can be found in the Supporting Information.

The O5′-H···G(N3) H-bond is not relevant for nonterminal nucleotides lacking a free 5′-OH terminal group. However, it could potentially be formed in some computations (especially MD simulations) where truncation of the system by removing the preceding DNA chain creates a 5′-terminal G with a 5′-OH terminal group. Because the initial distance between the O5′ and G(N3) atoms after such truncation is quite large, the O5′H···G(N3) H-bond was never formed in the geometry optimizations presented herein, which simplifies the analyses. The N2−H···O4′ internucleotide H-bond requires the first guanine to be in an anti conformation (Figure 3).44 This Hbond does not exist in real molecules because it is outcompeted by solvent molecules, and its formation in computations always

(PDB id:148D). The experimental structure contains one terminal syn G. However, the O5′-H···G(N3) H-bond is not unambiguously seen, and the deposited structure rather suggests that O5′ functions as the acceptor for the guanine amino group (i.e., this structure features an N2−H···O5′ interaction). While this conclusion may be due to an error in the refinement protocol, we decided not to polish the experimental data. In this particular case, the N2−H···O5′ Hbond was established in the QMrest and QMrest → QM/MM structures, whereas QM/MM and QM/MM → QM rest protocols located the more plausible O5′-H···G(N3) interaction. 2005

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation

Table 3. Differences (ΔΔE) between QM COSMO and MM PB ΔEs from Table 2 for the QMunrest, QM/MM, QMrest, QM/MM → QMrest, and QMrest → QM/MM Geometriesa bsc0

AA

AS

SA-aaaa

SA-aaab

SA-aabb

SA-abab

SA-abab-2

3AA+1SS

ΔΔEQMunrest ΔΔEQM/MM ΔΔEQMrest ΔΔEQM/MM→QMrest ΔΔEQMrest→QM/MM bsc0χOL4

0.0 0.0 0.0 0.0 0.0 AA

−4.9 0.0 −6.0 −10.0 −1.6 AS

+2.8 +22.5 +2.5 +5.0 +2.2 SA-aaaa

+7.4 +13.1 +8.3 +3.0 +9.4 SA-aaab

+2.8 +15.0 +4.6 0.0 +1.8 SA-aabb

+8.9 +1.0 +7.2 -2.9 +12.3 SA-abab

+4.2 +18.0 +5.4 +5.7 +16.3 SA-abab-2

−5.5 −0.2 +1.5 −4.3 +1.9 3AA+1SS

ΔΔEQMunrest ΔΔEQM/MM ΔΔEQMrest ΔΔEQM/MM→QMrest ΔΔEQMrest→QM/MM bsc0χOL4εζOL1

0.0 0.0 0.0 0.0 0.0 AA

−0.6 +3.0 −4.1 −5.8 +1.4 AS

+2.2 +24.5 +1.4 +5.8 +3.2 SA-aaaa

+9.7 +16.3 +10.3 +5.9 +12.2 SA-aaab

+4.6 +18.5 +5.0 +2.3 +3.8 SA-aabb

+8.7 +4.5 +7.9 −0.2 +13.6 SA-abab

+7.2 +21.0 +5.8 +7.5 +17.5 SA-abab-2

−12.6 −1.4 −2.0 −6.9 −1.4 3AA+1SS

ΔΔEQMunrest ΔΔEQM/MM ΔΔEQMrest ΔΔEQM/MM→QMrest ΔΔEQMrest→QM/MM

0.0 0.0 0.0 0.0 0.0

−0.4 +4.0 −1.8 −4.7 +2.9

+1.1 +24.3 +2.1 +6.1 +2.9

+8.4 +15.7 +10.6 +5.9 +11.8

+4.2 +17.1 +5.3 +2.6 +3.8

+7.4 +4.1 +9.7 −0.2 +13.8

+5.7 +20.9 +6.1 +7.6 +17.1

−11.1 −0.4 +0.9 −5.0 0.0

Calculated as ΔEQM − ΔEMM. Bold numbers indicate results that were used to compute the averaged correction for the force field’s inaccuracy when assessing the relative stability of the nonterminal syn-anti GpG step (see the main text for explanation).

a

obtained from single point DFT-D3 def2-TZVP COSMO calculations. The table also includes the corresponding MM data obtained from PB optimizations using three AMBER force field versions (see Methods) and starting from the QMunrest, QM/MM, QMrest, QM/MM → QMrest, and QMrest → QM/ MM geometries. Table 3 summarizes the differences between the QM and MM results. We first consider the QM data. There are several factors that do not affect the basic results. The relative energies computed for the QMrest K+ and QMrest Na+ geometries agree within 1 kcal/mol (cf. the fifth and last lines of Table 2). Changing the basis set from def2-TZVP to TZVP has only a marginal effect (∼3 kcal/mol on average). However, changing the continuum solvent model produces greater uncertainty. Replacing COSMO solvation with SMD (Table S1 in the Supporting Information) does not generally affect the basic trends, but the largest difference between the results for the two solvation models (8.4 kcal/mol, obtained for the QM/MM SA-aaaa structure optimized in the presence of K+) is quite significant. This demonstrates the substantial and unavoidable uncertainties introduced into energy calculations on nucleic acids by using continuum solvent models, which arise from the approximations these models use to compute solvation energies.134 More pronounced differences are observed between the relative energies calculated for the K+ and Na+ QM/MMoptimized structures. The most significant difference of 12 kcal/mol is found for the SA-aaaa system, for which the relative energies of the optimized K+ and Na+ structures are predicted to be +20 and +8 kcal/mol, respectively. The most striking differences emerge when comparing the relative energies obtained with QMrest and QM/MM optimizations. The most evident outlier is the SA-aaaa system: the relative energy with K+ channel ion is −24.9 kcal/mol when using the QMrest geometry and +20 kcal/mol when using the QM/MM geometry. Fortunately, the largest of these energy differences can be explained by considering the different numbers of O5′-H··· G(N3) terminal H-bonds (see above) in the corresponding

indicates severe non-native structural deformation. It can contribute erroneous additional stabilization and is associated with visible geometrical distortions, in particular unusually high β backbone angles, changes of helical twist, and quartet buckling. In previous work, these spurious N2−H···O4′ interactions were eliminated by introducing explicit water molecules, which also improved the β angles.44 Finally, weak C8−H···O5′ interactions may occur in purine nucleotides in certain geometries, as may C6−H···O5′ interactions in pyrimidines.131−133 These intramolecular interactions were not monitored in the preceding study.44 Such H-bonds are probably not very important in complete nucleic acid structures: although the relevant atoms are in moderately close proximity in some DNA X-ray structures, they are typically too widely separated to be considered H-bonded. Their frequent occurrence in earlier QM studies on nucleotides was probably a consequence of using small model systems. However, transient thermal sampling of such H-bonded configurations cannot be ruled out even in real structures (see below). Since these interactions can affect the computed data, we suggest that their presence/absence should be taken into account when analyzing the relative energies of structural isomers. Relative Energies of the Different Stems. In the following text, we first present the relative energies calculated on the basis of the QM/MM, QMrest, QM/MM → QMrest, and QMrest → QM/MM geometries and compare them to our previous44 estimates based on unrestrained QM optimized geometries (QMunrest). Then we compare the QM and MM energies. Finally, we comment on the main structural features of the optimized structures, i.e. their backbone angles, cation positions, cavity sizes, and Hoogsteen H-bonding patterns as well as the presence or absence of the H-bonding interactions discussed above. Table 2 presents the relative energies of the studied topologies with respect to the reference AA topology for each of the six series of optimized geometries, i.e. QMunrest (ref 44), QMunrest* with explicit waters (ref 44), QMrest, QM/MM, QM/ MM → QMrest, and QMrest → QM/MM. The energies are 2006

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation geometries. There are four such H-bonds in both the QMrest K+ and Na+ SA-aaaa structures. However, there is no such H-bond in the QM/MM K+ SA-aaaa structure and only one in the SAaaaa Na+ structure. Similarly, a single such H-bond occurs in the SA-aabb and SA-abab-2 QMrest structures but not in the equivalent QM/MM structures. The SA-abab QMrest structure features one presumably weaker (structurally less optimal) N2−H···O5′ interaction, while the QM/MM structure contains an O5′-H···G(N3) interaction. If we account for the estimated gain in potential energy of ∼5−9 kcal/mol per O5′-H···G(N3) H-bond,44 the major differences between the QM/MM and QMrest energies become dramatically smaller. An alternative table in which the calculated energies are modified by subtracting −7 kcal/mol from the relative energy of each structure for each such H-bond it contains can be found in the Supporting Information (Table S2). These H-bonds are absent in the QM/MM structures because the 5′-OH group hydrogen was oriented away from the N3 atom in the corresponding initial structures. The role of these H-bonds is clearly shown by the results of the second round of QM/MM optimizations, i.e. the QM/MM geometry reoptimizations starting from the QMrest geometries. These optimizations preserve all of the H-bonds established in the QMrest structures (see Supporting Information Table S10 for geometric details). The QMrest → QM/MM and QMrest structures thus contain the same sets of H-bonds, and their relative energy data are very similar (Table 2). In the opposite scenario, i.e. restrained optimization starting from the QM/MM-optimized geometries, only two terminal Hbonds out of a possible four are formed in the SA-aaaa structure. The relative energy for the resulting geometry is −11.3 kcal/mol, intermediate between those of the corresponding QMrest or QMrest → QM/MM (which both have four Hbonds) and QM/MM (no H-bond) geometries. As noted above, structures lacking the O5′-H···G(N3) H-bonds are entirely relevant, as both arrangements correspond to experimentally observable structures. However, it is important to account for the effects of O5′-H···G(N3) H-bonds consistently when computing relative energies. When this is done, all of the electronic structure calculations yield qualitatively consistent results. We stress that QM potential energy estimates such as those in Table 2 cannot be directly compared to experimental free energies. However, comparable free energy estimates can be obtained by combining the QM potential energy data with free energy estimates from explicit solvent MD simulations and using the QM calculations to correct the intrinsic inaccuracy of the MM force field. The procedure is outlined in ref 44. MM Computations. Table 3 presents differences between QM and MM relative energies for presumably equivalent geometries. The calculations are performed consistently at the QM//QM and MM//MM levels, i.e. the QM energies are computed on the QM-optimized geometries while the MM energies are computed on the MM-optimized geometries. All MM optimizations started from the corresponding QMoptimized geometries. The resulting data should be considered qualitatively rather than quantitatively. The individual energy points are affected by statistical (ensemble) uncertainties, which may explain some of the outliers. These uncertainties arise because we performed only one energy evaluation for each combination of stem topology, type of QM treatment, and force field dihedral parametrization variant: because the energy calculations can be

highly sensitive to the quality of the reference AA geometry, they will be strongly affected by the very limited sampling of the systems’ conformational space. To reduce this uncertainty, it would be necessary to repeat the calculations several times over, as is commonly done in (for example) the MM-PBSA method.57 An interesting but not unexpected result is that all three tested variants of the Cornell et al. force field give similar results. This indicates that conventional dihedral refinements of the Cornell et al. force field115 cannot alleviate the difference between QM and MM descriptions of the GQ stems. This is unsurprising because while adjustments of the dihedral terms significantly affect the behavior of nucleic acid simulations, they are not by themselves sufficient to achieve perfect tuning of force fields.2,43 Dihedral repametrizations may prevent the degradation of the simulated structures into certain clearly pathological geometries, such as untwisted canonical helices.2 However, they are not sufficient to eliminate the basic physical limitations of the pair-additive MM model.2 This suggestion is supported by our ongoing systematic benchmarking QM studies on smaller nucleic acid systems.9,10,38,55 Further tuning of the force field’s nonbonded terms could potentially improve the agreement between MM and QM results in future. However, it may be impossible to radically improve the agreement between QM and MM descriptions when using pairadditive MM methods with conformation-independent parameters such as constant point charges. In other words, the pair additive MM is likely to have principle accuracy limits that are not resolvable by reparametrizations. In this context, it should be noted that our earlier study demonstrated that the differences between the QM and MM results for the GQ stems were indeed due to the different descriptions of the DNA molecule rather than the solvent models used in each case.44 Another interesting result is that the MM reoptimizations eliminated the four spurious N2−H···O4′ H-bonds that heavily deformed the QMunrest AA and AS structures, thus effectively repairing these structures. This has several consequences. First, it makes the QMunrest and QMunrest → MMPB sets of geometries mutually structurally inconsistent. It should be noted that the ΔEQM − ΔEMM data for the QMunrest structures in Table 3 appear to be consistent with the QMrest data at first glance. However, this consistency could be entirely coincidental: the QMunrest and QMunrest → MMPB sets of structures are highly dissimilar and cannot be used as a basis for comparing QM and MM energies. Second, the improvement of the QMunrest → MMPB geometries might at first sight indicate that the MM description is more accurate than the QM description, because it avoids the formation of spurious H-bonds. This, however, would be an incorrect interpretation of the data. The MM model disrupts the spurious H-bonds because the MM description of the DNA structures is excessively rigid compared to the QM description. MM does not allow proper electronic structure redistributions, which contribute to stabilization of all H-bonds, and it does not provide enough flexibility around the β backbone torsion. This physical deficiency of the MM description in our particular case fortuitously compensates for the physical deficiency of the continuum treatment of the solvation. In contrast, the intrinsically more realistic QM description of the DNA molecule is not fully counterbalanced by the mean-field continuum solvent treatment and often develops gas-phase-like artifacts that would have to be suppressed by the inclusion of explicit water molecules. See also the discussion in ref 44 on page 9791. 2007

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation

We also note that the O6-cation distances (averaged over all O6-cation pairs) are systematically shorter in the QM/MM geometries than in the QMrest geometries for all topologies other than SA-aaab (Table S4 in the Supporting Information). Thus, the ions in the QM/MM structures are bound somewhat more tightly to the quartets than in the continuum solvent optimized structures. However, the difference is smaller than the standard deviations calculated for the eight O6-M+ distances. Similarly, the QM/MM → QMrest reoptimization leads to slightly looser binding of the channel ion. We also report data for the cavity sizes, measured in terms of the average distances between pairs of adjacent and opposite O6 atoms (Table S5 in the Supporting Information). The standard deviations presented in Tables S4 and S5 indicate that there is a greater degree of asymmetry in the QM/MM geometries across the strands of the stems, confirming that the QM/MM geometries are more heterogeneous than the QMrest ones. Backbone Angles−Restrained Optimization vs QM/MM Structures. The dihedral angles of the DNA backbone are important structural parameters, and specific DNA structures are partly defined by distinctive combinations of backbone dihedrals.9 Notably, many GQ topologies are intimately associated with noncanonical backbone substates. In the Supporting Information we document the backbone angles of all our starting and optimized structures (Figure S1 and Table S6). The restrained optimizations keep all backbone dihedrals close to their input values, as expected. The QM/MM approach also prevents the formation of non-native interactions that lead to high β angles, resulting in a remarkable improvement compared to the QMunrest geometries.44 However, specific starting configurations of water molecules can produce marked differences (i.e., heterogeneity) in the backbone angles of individual strands within a given GQ. This is particularly evident in the QM/MM-optimized geometries of the most symmetric topologies, AA and SA-aaaa, which exhibit pronounced variation in the β angles of their strands. In addition, unusual β angles are present in two strands of the QM/MM geometry of the 3AA+1SS topology with Na+ (a high β angle for strand 1 and a low β angle for strand 3). The AS QMrest → QM/MM geometry exhibits low β, γ, and δ+1 values in one strand. In addition, we found several cases in the QM/ MM geometries where some of the δ+1 angles are shifted to ∼80°, i.e. to the C3′-endo sugar pucker. As the C2′-endo and C3′-endo sugar puckers are close in energy, it is possible that a specific configuration of the water molecules can facilitate repuckering. Further investigations will be needed to better understand this behavior. Sugar puckering may be sensitive to the water environment or the truncation of the systems. Moreover, an accurate description of the sugar-pucker may be challenging even for QM methods. In principle, the uncertainty associated with structural heterogeneity of the QM/MM structures could be eliminated by performing a larger number of optimizations with randomized solvent configurations to mimic a more uniform (or averaged) solvent environment. It should be emphasized that the fact that the restrained optimizations remain closer to the starting structures than the noisier QM/MM computations does not mean that continuum solvent restrained optimization is universally better than QM/ MM optimization. The restrained optimizations are affected by uncertainties in the experimental structures, which can be sizable. Note, for example, the differences between different strands of those experimental structures that should ideally

Table 3 confirms a systematic difference between the QM and MM descriptions.44 Specifically, the SA stems with four syn-anti steps are typically destabilized by QM compared to MM. This table shows the ΔEQM − ΔEMM difference (ΔΔE) for 25 SA stem structures, each computed with three force field versions. These ΔEQM − ΔEMM values can be used to estimate the force field error of the earlier bsc0-based MM-PBSA calculation of the relative stability of the nonterminal syn-anti step with respect to the anti-anti step.57 However, some of the SA stems must be excluded from this computation. First, for reasons explained in the preceding paragraph, we have to exclude all data for the QMunrest set of geometries. Second, since the calculations primarily aim to derive the ΔEQM − ΔEMM difference for the common syn-anti step lacking the 5′-terminal H-bond, we also exclude stems with two or more such interactions, i.e., all SA-aaaa structures other than the QM/MM one. This is because there could be an additional nonequivalence between the MM and QM descriptions specifically associated with this terminal H-bond, which could bias the data. However, we tolerate SA stems with one terminal interaction (cf. Table 2) since otherwise we would have to exclude most of the structures. Finally, we use only data obtained with the bsc0 force field; to include the other force field variants rigorously, it would be necessary to repeat the MM-PBSA computations.57 (Nevertheless, the differences among different force fields versions appear to be small.) After accounting for these criteria, 17 SA stems were available for the computation; the data for these stems are highlighted in bold in Table 3. Their ΔEQM − ΔEMM differences range from −2.9 to 22.5 kcal/mol, which clearly shows the magnitude of the sampling noise in the computations. However, the average difference is 8.3 kcal/mol, which corresponds to 2.1 kcal/mol per GpG syn-anti step. This is the force field’s estimated overstabilization of a nonterminal syn-anti step relative to an anti-anti step. Our earlier MM-PBSA study estimated that the nonterminal syn-anti step is ∼3.5 kcal/ mol more stable than the anti-anti dinucleotide.57 Incorporating the 2.1 kcal/mol QM-MM correction would make the syn-anti and anti-anti dinucleotides essentially isoenergetic within the error margins of the used techniques, without any need for ad hoc computations to eliminate the spurious interactions.44 This is consistent with the empirical rule that GQs tend to maximize the number of anti-anti and syn-anti steps relative to the less stable anti-syn and syn-syn steps. Once again, we must stress that our calculations are based on very limited sampling. Note that although none of the considered structures contain any major structural artifacts, they are not homogeneous, and especially those obtained by the QM/MM protocol contain subtle deformations such as some quartet nonplanarities (buckling). Channel Ion Binding in Restrained and QM/MM Optimizations. We performed several additional analyses that are fully documented in the Supporting Information. Among other things, we re-evaluated the relative energies for the QM/ MM and QMrest geometries using single point COSMO DFTD3 energy calculations without the channel ions (Table S3 in the Supporting Information). Removing the ions caused only minor changes (∼3 and ∼4 kcal/mol for QM/MM and QMrest geometries, respectively) in the relative stabilities compared to Table 2. This confirms two things. First, the minor heterogeneities in ion positions that occur in some computations do not affect the basic results. Second, the relative energy variability reported in our study is dominated by the internal DNA conformational energy. 2008

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation have identical geometries, such as the AA stem. Non-negligible bias can be introduced into the starting structures by the inherently less accurate NMR data. A textbook example is the 2GKU102 NMR structure, which contains some G-nucleotides with a β angle of ∼300°. Such β angle values are highly unusual (not known from X-ray structures) and may be a result of errors introduced by the NMR structure refinement protocol. One such high-β angle is transferred to our computations (in strand 3 of the 3AA+1SS structure). Of course, in principle, we could leave one angle unrestrained or change its target value. However, due to the genuine correlation between the backbone angles, we do not think this would be a generally applicable procedure. Perhaps, when working from a high-resolution experimental structure and using just a single optimization per structure, restrained optimization would be a more conservative approach. However, as and when adequate computer resources become available, we would suggest performing an extended series of QM/MM optimizations with diverse initial configurations to better map the native region of the conformational space. The QM/MM protocol is in principle suitable for obtaining an ensemble of structures that could enable better characterization of the potential energy surfaces of the DNA building blocks. Nevertheless, both restrained and QM/MM optimizations appear to efficiently prevent the formation of spurious structures such as those biased by the N2−H···O4′ Hbonds, and the two approaches complement each other. Intramolecular H-Bonds. Let us make a few comments on the H-bonds. All methods and topologies preserve Hoogsteen H-bond pairings. We do not observe any large deformations or bifurcated arrangements of these interactions. However, there are some notable minor features: (i) larger variance is observed for the QM/MM geometries, as already discussed for the backbone angles; and (ii) the N7···H bonds are systematically shorter in structures featuring Na+ compared to those with K+. The H-bond lengths of the Hoogsteen pairings are summarized in Table S7 of the Supporting Information. The C8−H···O5′ interactions (Figure 3) are weak, but a twoquartet AA stem can potentially form 8 such interactions. The O5′···C8 and the respective O5′···H8 distances for the AA, 3AA+1SS, and AS topologies are presented in Table S8 of the Supporting Information. The data indicate an absence of such interactions in the parent experimental structures. In contrast, many computed geometries indicate possible weakly stabilizing O5′···H8 contacts at the 5′-termini in all three topologies. It is interesting to note that the QMunrest* AA geometry that was originally optimized with explicit waters to suppress the N2− H···O4′ interaction is the only one to have short O5′···H8 distances at both the 3′- and 5′-termini. Thus, we cannot fully rule out the possibility that this contact is sampled to some extent. We also performed an explicit solvent classical MD simulation of the AA stem (Figure 4); a small fraction of the simulated thermal fluctuations produced O5′···C8 distances shorter than 3.5 Å. Details on this MD simulation are given in the Supporting Information section titled “Protocol for the MD simulation on the AA-K+ topology”. The N2−H···O4′ interaction is, as explained above, nonnative. This spurious interaction occurs only in the original QMunrest optimized geometries without explicit solvent molecules (Table S9 in the Supporting Information). Both QMrest and QM/MM procedures prevent this interaction. An exception is noted for strand 1 of the AS QM/MM geometry with Na+, though the H-bond is longer than in the QMunrest optimization. This indicates that, due to lack of restraints, the

Figure 4. Distribution of C8···O5′ distances for strands 1−4 (from left to right) in the AA structure in the course of a 2 ns explicit solvent MD simulation. Top: 5′-end. Bottom: 3′-end.

QM/MM procedure can occasionally form non-native interactions when certain initial solvent configurations are used. Therefore, inspection of the optimized QM/MM structures is always advised. For the sake of completeness, O5′-H···G(N3) (or, sometimes N2−H···O5′ in SA-abab) interactions in syn 5′-termini are observed in one strand of the SA-abab K+ and SA-aaaa Na+ topologies in the QM/MM geometries. In the QMunrest and QMrest geometry sets, this interaction forms in one strand of the SA-aabb, SA-abab, and SA-abab-2 topologies and in all four SAaaaa strands; in all such cases, it is preserved in the corresponding QMrest → QM/MM geometries. This approach thus yields QM/MM geometries featuring the interaction in all 5′-terminal syn nucleotides, with short interheteroatom distances. In the QM/MM → QMrest geometry sets, this interaction forms in two strands of the SA-aaaa topology as well as one strand of the SA-aabb, SA-abab, and SA-abab-2 topologies. These results are explained in detail above (cf. also their listings in Table 2). Additional Calculations on the AA Model. In the following paragraphs, we present a series of QM/MM calculations on the AA structure with K+ exploring the effects of various methodological parameters on the resulting structures. Specifically, we examine the potential bias of the starting geometry, the basis set of the QM region, the water model of the MM level, the number of solvation shells, the initial configuration of waters, the presence of counterions at the MM level, the application of electrostatic embedding (EE), and the potential benefits of freezing the outer layer of solvent molecules. We make only brief comments on selected results because none of the tests identified any specific weakness of our basic QM/MM protocol. The calculations are performed with different input parameters (calculations A1−A4, B1, B2, C1−C3, and D1, D2), the details of which are summarized in Table 4. The superscripted label ee denotes a calculation performed using EE. All “A” calculations have the four “bulk” counterions treated at the MM level. Calculation A1 is the same QM/MM calculation on the AA topology as in the first part of the article, and we use it as the reference calculation. The rest of the “A” calculations explore the effects of using EE (A1ee), a different starting geometry (A2) or basis set (A3), and the SPC/E135 water parametrization (A4, A4ee). The “B” calculations exclude the bulk counterions. The number of water molecules is determined by a short QM/MM MD simulation with a frozen solute molecule (for further details see the Supporting Information section on “Estimating 2009

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation

Table 4. Input Details of the Test QM/MM Calculation Sets A, B, C, and D, Which Were Performed on the AA GQ Stem models

stema

A1 A1ee A2 A3 A4 A4ee B1f B1eef B2f B3f B3eef C1 C1ee D1 D2

QM QM QM QM QM QM QM QM QM QM QM QM QM QM QM

waterb MM MM MM MM MM MM MM MM MM MM MM MM MM MM MM

(604) (604) (478) (604) (604) (604) (137) (137) (266) (416) (416) (604) (604) (964) (964)

counter ionsc

water model

basis set

EEd

input geom

MM (4) MM (4) MM (4) MM (4) MM (4) MM (4) - (0) - (0) - (0) - (0) - (0) QM (4) QM (4) QM (4) QM (4)

TIP3P TIP3P TIP3P TIP3P SPC/E SPC/E TIP3P TIP3P TIP3P TIP3P TIP3P TIP3P TIP3P TIP3P TIP3P

TZVP TZVP TZVP def2e TZVP TZVP TZVP TZVP TZVP TZVP TZVP TZVP TZVP TZVP TZVP

no yes no no no yes no yes no no yes no yes no no

X-ray A1 QMunrest44 A1 X-ray A4 X-ray B1 X-ray X-ray B3 X-ray C1 X-ray X-ray

a Method employed for the stem including the channel ion. bMethod employed for water molecules and their number (in parentheses). cMethod employed for the bulk counterions and number of them in parentheses. dElectrostatic embedding. edef2-TZVP. fThe number of waters included in models B1−B3 was estimated after short MD simulation (see paragraph entitled “Estimate of Solvation Shells by MD” in the Supporting Information for details).

Figure 5. Comparison of X-ray99 and QMunrest*44 backbone dihedral angles to those of models A1, B1, B2, and B3.

non-native interactions due to the inclusion of explicit solvent molecules at the MM level. The A2 QM/MM optimization started from the AA QMunrest-optimized structure with the spurious H-bonds. When the spurious H-bonds are present in the starting structure, the QM/MM procedure is not able to disrupt them. As discussed above, the QM/MM optimization is also not capable of forming the 5′-end O5′-H···G(N3) Hbonds if the hydroxyl hydrogen is initially not oriented toward the expected H-bond acceptor (although such structures are not invalid per se). Therefore, when performing QM/MM optimizations, one must always verify that the starting structures contain all of the desired H-bonds with no undesired interactions. The MM waters and gradient optimization can prevent the formation of H-bonds but not disrupt them if they are already present. Additional insights into the dependence of the QM/MM calculations on the starting structure can be obtained from the B1, B2, and B3 optimizations. These optimizations all lack bulk counterions but differ in the number of included waters and the starting water configuration. The results are rather similar to those for the reference calculation A1. The β angles are

Solvation Shells by MD” and Table S11); calculations are performed with one (B1, B1ee), two (B2), and three (B3, B3ee) solvation shells around the stem. The starting configurations of the water molecules differ in each case because the snapshots are taken from different instances of the simulation. Calculations C1 and C1ee include the counterions at the QM level, while all other calculation parameters are the same as for A1. Finally, calculation D1 has a large solvent shell (extending around 15 Å from the solute), of which the outer layer (which is approximately 5 Å thick) is kept frozen during the optimization. All of the input parameters of D2 are identical to those of D1 but without the outer layer being frozen. As mentioned in the Introduction, previous QM unrest geometry optimizations44 resulted in four spurious N2−H··· O4′ interactions with high β values (around 240°) in the AA stem. Addition of X-ray water molecules at suitable positions eliminated these interactions and brought the β angles closer to their experimental values although they remained overly high at around 200°. Our QM/MM optimization A1 starting from the same input geometry as in ref 44 successfully eliminates the 2010

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Journal of Chemical Theory and Computation

■

somewhat overestimated, being similar to those obtained in the QMunrest* optimization with four explicit water molecules, but never as large as in the QMunrest optimization with no waters (Figure 5 and Table S12 in the Supporting Information.) We do not see any systematic trends in these calculations, and all of them appear to shield the system from forming spurious interactions. As the optimizations differ with respect to the number and initial positions of the waters, the observed differences primarily reflect random differences in the starting solvent configurations. This indicates that even a reduced water shell provides a quite satisfactory environment for the present system. The A3 optimization was a repetition of the A1 calculation with the larger def2-TZVP basis set. Because of the increased computational cost due to the larger basis set, it was started from the optimized A1 geometry. The backbone angles of the A3 geometry were practically identical to those of the A1. The geometries obtained in the ME calculations using the TIP3P (A1) and SPC/E (A4) water models were not noticeably different. The same was true for the calculations with (D1) and without (D2) a frozen external layer of solvent. The most important part of our test computations was our comparison of mechanical and electrostatic embedding. Because the EE optimizations converged much more slowly than their ME counterparts, each EE optimization (A1ee, A4ee, B1ee, B3ee, and C1ee) started from the corresponding MEoptimized structure. The MUEs of the backbone angles with respect to the experimental structures (Table S13 in the Supporting Information) for each of the EE calculations were practically identical, although there was a slight tendency toward lower β angles compared to the ME-optimized structures in the B1ee and B3ee runs. Overall, though, none of the five EE calculations yielded results that differed visibly from the corresponding ME computations. However, it is often harder for EE geometry optimizations to converge, and the steps required to achieve convergence may differ from case to case. Given the lack of significant differences in the backbone dihedral angles among the different QM/MM geometry optimizations, we did not extend our analysis to other geometric features. However, for completeness we include data on cavity size, channel ion positions, and intramolecular Hbonds for this set of calculations in the Supporting Information Tables S14−S18. Systematic Total Electronic Energy Difference between QM/MM and QMrest Optimized Structures. We have noticed that, when used in single point QM COSMO energy calculations, all QM/MM geometries obtained in this study are higher in energy than the corresponding QMrest optimized geometries. We suggest that this is a cumulative effect of minor structural differences between the QM/MM and QMrest sets of structures, with the QM/MM structures being slightly less structurally relaxed than the QM COSMO structures. Further details on this topic are given in the last section of the Supporting Information (titled “Relative energies of QM and QM/MM geometries”). Although this does not introduce any error into the computations, it means that we should not mix QM/MM and QM-COSMO optimized structures in energy comparisons. The difference can be neatly eliminated if the QM/MM-optimized structure is subjected to a restrained COSMO QM optimization.

Article

CONCLUDING REMARKS

MD simulations are the most widely used computational tools for studying nucleic acids despite the known limitations of the MM force fields.1,2,136 Applications of more accurate QM methods have been limited to small model systems but have provided unique insights into the H-bonding, stacking, ion binding, and backbone conformations of nucleic acids. In recent years, advances in computer hardware and sophisticated theoretical approximations have significantly increased the capabilities and potential of both MD and QM methods. MD simulations have become increasingly long, nowadays routinely reaching the μs time scale for medium-sized nucleic acid structures.1,2,136 Moreover, the development of highly accurate dispersion-corrected DFT methods has for the first time made it possible to study small but complete nucleic acid building blocks using methods capable of reliably describing their electronic structure. QM studies can thus be extended to complete nucleic acid building blocks such as the two-quartet GQ stems studied here. However, because such studies remain limited to the description of potential energy surfaces, it will be important to find synergies with the MM-based approaches to obtain biochemically relevant data.44 The nucleic acid systems that are currently amenable to analysis using QM computations are not fully complete due to the need to use implicit solvent environments and to truncate the systems to a computationally affordable size. These issues were exemplified in our preceding work where several of the optimized two-quartet GQ stem structures were spoiled by non-native intramolecular H-bonds and spurious backbone conformations.44 Here, as an extension of our previous studies, we performed both restrained QM optimizations59 in continuum solvent and QM/MM optimizations in explicit solvent on a series of GQ stem topologies. Both approaches represent a clear step forward from unrestrained implicit solvent QM optimizations even though they have distinct weaknesses. Restrained optimizations reliably produce geometries close to the target geometries. However, if accurate target structures are not available, all of the uncertainties and errors in the target geometries (which may arise from the analysis of low-resolution X-ray and NMR structures) are inevitably transferred to the optimized geometries. In some cases, these errors could probably be reduced by knowledge-based modifications of the target structures. However, due to the known interdependency of individual nucleic acids’ backbone dihedrals,38 such rational interventions are not a panacea. On the other hand, while the inclusion of explicit solvent molecules by QM/MM is a realistic alternative, it is still too expensive to account for the dynamic evolution of the solvent. This makes it necessary to use a single specific solvent configuration in each optimization, which creates some heterogeneity in the resulting geometries, as discussed above. We suggest that further regularization of the data can be achieved by successively employing QM/MM and restrained QM optimizations. Both restrained QM continuum solvent and QM/MM explicit solvent optimizations prevent the formation of nonnative H-bonds, namely the N2−H···O4′ guanosine interactions. We nevertheless suggest that H-bonding should always be carefully monitored in the calculations. The QM/MM protocol may occasionally fail to prevent the formation of a spurious H-bond due to a specific initial water configuration, or it may hinder the formation of a desired intramolecular H-bond 2011

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation

applied to structures sampled during MD simulations to achieve a refined analysis of the potential energy surfaces of nucleic acid building blocks. Such calculations can be used to better interpret experimental data, to compensate for errors introduced by MM force fields, and to verify and parametrize other computational methods. We suggest that continuing improvements in the capabilities of hardware and software, together with the ongoing development of fast QM methods, will in the near future make it possible to routinely perform a series of hundreds to thousands of individual QM computations on biomolecular building blocks with far greater accuracy than is possible with an MM description. Such computations, when applied with smart strategies for conformational space sampling and combined with additional information from MM methods, should substantially increase our ability to derive robust theoretical descriptions of the chemistry and structures of nucleic acids. It is however important to note that for any biochemically relevant application of QM methods to nucleic acids, we need to ensure that the calculations are performed on realistic geometries. The use of deformed geometries would wipe out any advantage of the higher intrinsic accuracy of the QM methods compared to the MM methods. Obviously, in an ideal case, computational chemistry should investigate the studied systems using computations that are not biased by experimental data. Unfortunately, biomolecular systems are too complex and delicate to carry out computations that are entirely independent of experimental data. Contemporary computational chemistry is unable (and will not become able for the foreseeable future) to predict biomolecular structures and energies from unbiased computations. Careful utilization of the available experimental structural data is especially important in QM studies due to their very limited sampling compared to MM-based methods.

due to the initial orientation of the relevant hydrogen. It is to be noted that various kinds of restrained optimizations are commonly used in connection with MM approaches for several decades. However, to the best of our knowledge, this approach is not commonly used in QM studies of biomolecular systems and is not standardly available in many QM codes. We suggest that the use of appropriate soft structural restraints is imperative in QM studies of biomolecular building blocks, as unconstrained QM optimizations typically result in structures that lack biochemical relevance and are unsuitable for further analyses. Standard hard constraints are too severe, and their application often results in “off-target” deformations and energy strains. It should also be pointed out that in order to get energies using single-point continuum solvent QM calculations after QM/MM geometry optimizations, we remove the explicit water molecules. This may affect the resulting energies by eliminating some specific water contacts. Because we are primarily interested in comparing the relative energies of different stem conformers, this approach seems justifiable. However, it is always necessary to monitor the presence and absence of individual intramolecular H-bonds in the studied structures, as is done in the Results and Discussion section of this paper. If two structures have different sets of H-bonds, this must be accounted for when comparing and interpreting their relative energies. Our improved results qualitatively confirm the previously reported findings based on QM calculations of the relative energies of different GQ stem topologies using unrestrained QM optimizations in implicit solvent environments44 and greatly reduce the associated uncertainty. We also explored the dependence of the results on various QM/MM methodological parameters. We tested the effect of the starting geometry, the MM water model, the basis set of the QM region, the number of water molecules included in the model and their initial configuration, the inclusion of counterions at the MM or QM levels, the effect of freezing an outer solvent layer, and most importantly the effect of electrostatic embedding. None of the tested factors affected the results. Let us specifically comment on the use of mechanical embedding. The electrostatic embedding is usually considered as physically more realistic than the mechanical embedding in QM/MM calculations. However, as shown above, for our particular purpose both mechanical and electrostatic embedding lead to essentially identical results. Although this might be at first sight surprising, there are several reasons that explain this observation. In our approach, we immerse a chemically complete biomolecular building block that is described by QM into a box of explicit MM waters. Thus, the QM/MM border does not cut through any covalent bonds, and there are no link atoms involved in our QM/MM computations. Further, the MM environment is essentially uniform. Aforementioned factors minimize inhomogeneous polarization effects that otherwise would require electrostatic embedding and justify the use of mechanical embedding in our work, as directly confirmed by a series (five structures) of our comparative computations. Note also that the energies are not calculated using the QM/MM protocol. In summary, both restrained continuum solvent QM and QM/MM explicit solvent optimizations provide significantly more stable results than standard QM optimizations in structural relaxations of nucleic acid building blocks. The two methods are complementary and should perform well when

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b01025. List of abbreviations; relative energies of all topologies using the SMD method; relative energies of all topologies adjusted for the O5′-H···G(N3) H-bonds; relative energies of all topologies without channel ions; O6···K+ distances; distances between adjacent and opposite O6 atoms; schematic representation of backbone dihedral angle profiles; list of backbone dihedral angles; N7···H21 and O6···H1 H-bond distances; O5′···C8 and O5′···H8 distances; N2(G)···O4′ (G+1) and H22(G)···O4′ (G+1) distances; O5′···N3 and H···N3 distances; summary of MD protocol on AA-K+; complete structural and energy data on the additional set of QM/MM calculations on the AA topology (PDF) Cartesian coordinates of optimized geometries of all topologies (ZIP)

■

AUTHOR INFORMATION

Corresponding Author

*Phone: 420 541 517 133. Fax: 420 541 211 293. E-mail: [email protected]. Notes

The authors declare no competing financial interest. 2012

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation

■

(15) Gu, J.; Leszczynski, J. Origin of Na/K selectivity of the guanine tetraplexes in water: The theoretical rationale. J. Phys. Chem. A 2002, 106, 529−532. (16) Louit, G.; Hocquet, A.; Ghomi, M.; Meyer, M.; Sühnel, J. Are guanine tetrads stabilised by bifurcated hydrogen bonds? An AIM topological analysis of the electronic density. PhysChemComm 2002, 5, 94−98. (17) Meyer, M.; Steinke, T.; Brandl, M.; Sühnel, J. Density functional study of guanine and uracil quartets and of guanine quartet/metal ion complexes. J. Comput. Chem. 2001, 22, 109−124. (18) van Mourik, T.; Dingley, A. J. Characterization of the monovalent ion position and hydrogen-bond network in guanine quartets by DFT calculations of NMR parameters. Chem. - Eur. J. 2005, 11, 6064−6079. (19) Gu, J.; Leszczynski, J. A remarkable alteration in the bonding pattern: An HF and DFT study of the interactions between the metal cations and the Hoogsteen hydrogen-bonded G-tetrad. J. Phys. Chem. A 2000, 104, 6308−6313. (20) Louit, G.; Hocquet, A.; Ghomi, M.; Meyer, M.; Sühnel, J. Guanine tetrads interacting with metal ions. An AIM topological analysis of the electronic density. PhysChemComm 2003, 6, 1−5. (21) Gresh, N.; Pullman, B. A theoretical study of the selective entrapment of alkali and ammonium cations between guanine tetramers. Int. J. Quantum Chem. 1985, 28, 49−56. (22) Yurenko, Y. P.; Novotný, J.; Sklenár,̌ V.; Marek, R. Exploring Non-Covalent Interactions in Guanine- and Xanthine-Based Model DNA Quadruplex Structures: a Comprehensive Quantum Chemical Approach. Phys. Chem. Chem. Phys. 2014, 16, 2072−2084. (23) 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. (24) Chawla, M.; Abdel-Azeim, S.; Oliva, R.; Cavallo, L. Higher order structural effects stabilizing the reverse watson-crick guanine-cytosine base pair in functional RNAs. Nucleic Acids Res. 2014, 42, 714−726. (25) Wilson, K. A.; Kellie, J. L.; Wetmore, S. D. DNA-protein piinteractions in nature: Abundance, structure, composition and strength of contacts between aromatic amino acids and DNA nucleobases or deoxyribose sugar. Nucleic Acids Res. 2014, 42, 6726−6741. (26) Halder, S.; Bhattacharyya, D. RNA structure and dynamics: A base pairing perspective. Prog. Biophys. Mol. Biol. 2013, 113, 264−283. (27) Wells, R. A.; Kellie, J. L.; Wetmore, S. D. Significant strength of charged DNA-protein π-π Interactions: A preliminary study of cytosine. J. Phys. Chem. B 2013, 117, 10462−10474. (28) Lech, C. J.; Phan, A. T.; Michel-Beyerle, M. E. Electron-hole transfer in G-quadruplexes with different tetrad stacking geometries: A combined QM and MD study. J. Phys. Chem. B 2013, 117, 9851−9856. (29) 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. (30) Hongo, K.; Cuong, N. T.; Maezono, R. The importance of electron correlation on stacking interaction of adenine-thymine basepair step in B-DNA: A quantum Monte Carlo study. J. Chem. Theory Comput. 2013, 9, 1081−1086. (31) Hill, J. G.; Platts, J. A. Calculating stacking interactions in nucleic acid base-pair steps using spin-component scaling and local second order M?ller-Plesset perturbation theory. Phys. Chem. Chem. Phys. 2008, 10, 2785−2791. (32) Cooper, V. R.; Thonhauser, T.; Puzder, A.; Schröder, E.; Lundqvist, B. I.; Langreth, D. C. Stacking interactions and the twist of DNA. J. Am. Chem. Soc. 2008, 130, 1304−1308. (33) Leininger, M.; Nielsen, I. M. B.; Colvin, M. E.; Janssen, C. L. Accurate structures and binding energies for stacked uracil dimers. J. Phys. Chem. A 2002, 106, 3850−3854. (34) Luo, R.; Gilson, H. S. R.; Potter, M. J.; Gilson, M. K. The physical basis of nucleic acid base stacking in water. Biophys. J. 2001, 80, 140−148. (35) Fonseca Guerra, C.; Zijlstra, H.; Paragi, G.; Bickelhaupt, M. Telomere Structure and Stability: Covalency in Hydrogen Bonds, Not

ACKNOWLEDGMENTS This work was funded by the Czech Science Foundation (P208/11/1822) and by the Ministry of Education, Youth and Sports of the Czech Republic under the project CEITEC 2020 (LQ1601). We also acknowledge support by Praemium Academiae from the Academy of Sciences of the Czech Republic. We would like to thank one of our Reviewers for exceptionally careful reading of the paper and many highly relevant comments.

■

REFERENCES

(1) Cheatham, T. E., III; Case, D. A. Twenty-five years of nucleic acid simulations. Biopolymers 2013, 99, 969−977. (2) Šponer, J.; Banás,̌ P.; Jurečka, P.; Zgarbová, M.; Kührová, P.; Havrila, M.; Krepl, M.; Stadlbauer, P.; Otyepka, M. Molecular dynamics simulations of nucleic acids. From tetranucleotides to the ribosome. J. Phys. Chem. Lett. 2014, 5, 1771−1782. (3) Burda, J. V.; Sponer, J.; Leszczynski, J.; Hobza, P. Interaction of DNA Base Pairs with Various Metal Cations (Mg2+, Ca2+, Sr2+, Ba2+, Cu+, Ag+, Au+, Zn2+, Cd2+, And Hg2+): Nonempirical Ab Initio Calculations on Structures, Energies, and Nonadditivity of the Interaction. J. Phys. Chem. B 1997, 101, 9670−9677. (4) Hobza, P.; Šponer, J. Structure, Energetics, and Dynamics of the Nucleic Acid Base Pairs: Nonempirical Ab Initio Calculations. Chem. Rev. 1999, 99, 3247−3276. (5) Hobza, P.; Šponer, J. Toward True DNA Base-Stacking Energies: MP2, CCSD(T), and Complete Basis Set Calculations. J. Am. Chem. Soc. 2002, 124, 11802−11808. (6) Mignon, P.; Loverix, S.; Steyaert, J.; Geerlings, P. Influence of the Pi-Pi Interaction on the Hydrogen Bonding Capacity of Stacked DNA/RNA Bases. Nucleic Acids Res. 2005, 33, 1779−1789. (7) Šponer, J.; Jurečka, P.; Marchan, I.; Luque, F. J.; Orozco, M.; Hobza, P. Nature of Base Stacking: Reference Quantum-Chemical Stacking Energies in Ten Unique B-DNA Base-Pair Steps. Chem. - Eur. J. 2006, 12, 2854−2865. (8) Meyer, M.; Hocquet, A.; Sühnel, J. Interaction of Sodium and Potassium Ions with Sandwiched Cytosine-, Guanine-, Thymine-, and Uracil-Base Tetrads. J. Comput. Chem. 2005, 26, 352−364. (9) Mládek, A.; Krepl, M.; Svozil, D.; Cech, P.; Otyepka, M.; Banás,̌ P.; Zgarbová, M.; Jurečka, P.; Šponer, J. Benchmark QuantumChemical Calculations on a Complete Set of Rotameric Families of the DNA Sugar−Phosphate Backbone and Their Comparison with Modern Density Functional Theory. Phys. Chem. Chem. Phys. 2013, 15, 7295−7310. (10) Mládek, A.; Banás,̌ P.; Jurečka, P.; Otyepka, M.; Zgarbová, M.; Šponer, J. Energies and 2′-Hydroxyl Group Orientations of RNA Backbone Conformations. Benchmark CCSD(T)/CBS Database, Electronic Analysis, and Assessment of DFT Methods and MD Simulations. J. Chem. Theory Comput. 2014, 10, 463−480. (11) Gkionis, K.; Kruse, H.; Platts, J. A.; Mládek, A.; Koča, J.; Šponer, 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. (12) Lech, C. J.; Heddi, B.; Phan, A. T. Guanine base stacking in Gquadruplex nucleic acids. Nucleic Acids Res. 2013, 41, 2034−2046. (13) Fonseca Guerra, C.; van der Wijst, T.; Poater, J.; Swart, M.; Bickelhaupt, F. M. Adenine versus guanine quartets in aqueous solution: Dispersion-corrected DFT study on the differences in pstacking and hydrogen-bonding behavior. Theor. Chem. Acc. 2010, 125, 245−252. (14) Otero, R.; Schöck, M.; Molina, L. M.; Lægsgaard, E.; Stensgaard, I.; Hammer, B.; Besenbacher, F. Guanine quartet networks stabilized by cooperative hydrogen bonds. Angew. Chem., Int. Ed. 2005, 44, 2270−2275. 2013

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation Resonance Assistance, Causes Cooperativity in Guanine Quartets. Chem. - Eur. J. 2011, 17, 12612−12622. (36) Yurenko, Y. P.; Novotný, J.; Mitoraj, M. P.; Sklenár,̌ 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. (37) Foloppe, N.; MacKerell, A. D., Jr. Contribution of the phosphodiester backbone and glycosyl linkage intrinsic torsional energetics to DNA structure and dynamics. J. Phys. Chem. B 1999, 103, 10955−10964. (38) Šponer, J.; Mládek, A.; Šponer, J. E.; Svozil, D.; Zgarbová, M.; Banás,̌ P.; Jurečka, P.; Otyepka, M. The DNA and RNA sugarphosphate backbone emerges as the key player. An overview of quantum-chemical, structural biology and simulation studies. Phys. Chem. Chem. Phys. 2012, 14, 15257−15277. (39) Zgarbová, M.; Luque, F. J.; Šponer, J.; Cheatham, T. E., III; Otyepka, M.; Jurečka, P. Toward improved description of DNA backbone: Revisiting epsilon and zeta torsion force field parameters. J. Chem. Theory Comput. 2013, 9, 2339−2354. (40) Perez, A.; Marchan, I.; Svozil, D.; Šponer, J.; Cheatham, T. E., III; Laughton, C. A.; Orozco, M. Refinement of the AMBER Force Field for Nucleic Acids: Improving the Description of α/γ Conformers. Biophys. J. 2007, 92, 3817−3829. (41) Krepl, M.; Zgarbova, M.; Stadlbauer, P.; Otyepka, M.; Banas, P.; Koca, J.; Cheatham, T. E., III; Jurecka, P.; Sponer, J. Reference Simulations of Noncanonical Nucleic Acids with Different χ Variants of the AMBER Force Field: Quadruplex DNA, Quadruplex RNA, and Z-DNA. J. Chem. Theory Comput. 2012, 8, 2506−2520. (42) Hart, K.; Foloppe, N.; Baker, C. M.; Denning, E. J.; Nilsson, L.; MacKerell, A. D., Jr. Optimization of the CHARMM Additive Force Field for DNA: Improved Treatment of the BI/BII Conformational Equilibrium. J. Chem. Theory Comput. 2012, 8, 348−362. (43) Šponer, J.; Cang, X.; Cheatham, T. E., III Molecular Dynamics Simulations of G-DNA and Perspectives on the Simulation of Nucleic Acid Structures. Methods 2012, 57, 25−39. (44) Šponer, J.; Mládek, A.; Špačková, N.; Cang, X.; Cheatham, T. E., III; 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. (45) Zubatiuk, T. A.; Shishkin, O. V.; Gorb, L.; Hovorun, D. M.; Leszczynski, J. B-DNA Characteristics Are Preserved in Double Stranded D(A)3·D(T)3 and D(G)3·D(C)3 mini-Helixes: Conclusions From DFT/M06-2X Study. Phys. Chem. Chem. Phys. 2013, 15, 18155− 18166. (46) Kruse, H.; Havrila, M.; Šponer, J. QM computations on complete nucleic acids building blocks: Analysis of the Sarcin-Ricin RNA motif using DFT-D3, HF-3c, PM6-D3H, and MM approaches. J. Chem. Theory Comput. 2014, 10, 2615−2629. (47) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456−1465. (48) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (49) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215−241. (50) Zhao, Y.; Truhlar, D. G. Applications and validations of the Minnesota density functionals. Chem. Phys. Lett. 2011, 502, 1−13. (51) Grimme, S. Density functional theory with London dispersion corrections. Wiley Interdiscip. Rev. Mol. Sci. 2011, 1, 211−228. (52) Goerigk, L.; Grimme, S. Efficient and Accurate Double-HybridMeta-GGA Density FunctionalsEvaluation with the Extended GMTKN30 Database for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2011, 7, 291−309.

(53) Goerigk, L. How Do DFT-DCP, DFT-NL, and DFT-D3 Compare for the Description of London-Dispersion Effects in Conformers and General Thermochemistry? J. Chem. Theory Comput. 2014, 10, 968−980. (54) Goerigk, L.; Kruse, H.; Grimme, S. Benchmarking Density Functional Methods against the S66 and S66 × 8 Datasets for NonCovalent Interactions. ChemPhysChem 2011, 12, 3421−3433. (55) Kruse, H.; Mladek, A.; Gkionis, K.; Hansen, A.; Grimme, S.; Sponer, J. Quantum Chemical Benchmark Study on 46 RNA Backbone Families Using a Dinucleotide Unit. J. Chem. Theory Comput. 2015, 11, 4972−4991. (56) Peverati, R.; Truhlar, D. G. Quest for a universal density functional: the accuracy of density functionals across a broad spectrum of databases in chemistry and physics. Philos. Trans. R. Soc., A 2014, 372, 20120476. (57) Cang, X.; Sponer, J.; Cheatham, T. E., III Explaining the Varied Glycosidic Conformational, G-Tract Length and Sequence Preferences for Anti-Parallel G-Quadruplexes. Nucleic Acids Res. 2011, 39, 4499− 4512. (58) Sponer, J.; Šponer, J. E.; Mládek, A.; Banás,̌ P.; Jurečka, P.; Otyepka, M. How to Understand Quantum Chemical Computations on DNA and RNA Systems? a Practical Guide for Non-Specialists. Methods 2013, 64, 3−11. (59) Kruse, H.; Sponer, J. Towards biochemically relevant QM computations on nucleic acids: controlled electronic structure geometry optimization of nucleic acid structural motifs using penalty restraint functions. Phys. Chem. Chem. Phys. 2015, 17, 1399−1410. (60) Burge, S.; Parkinson, G. N.; Hazel, P.; Todd, A. K.; Neidle, S. Quadruplex DNA: sequence, topology and structure. Nucleic Acids Res. 2006, 34, 5402−5415. (61) Neidle, S. The structures of quadruplex nucleic acids and their drug complexes. Curr. Opin. Struct. Biol. 2009, 19, 239−250. (62) Li, J.; Correia, J. J.; Wang, L.; Trent, J. O.; Chaires, J. B. Not so crystal clear: the structure of the human telomere G-quadruplex in solution differs from that present in a crystal. Nucleic Acids Res. 2005, 33, 4649−4659. (63) 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. (64) Adrian, M.; Heddi, B.; Phan, A. T. NMR spectroscopy of Gquadruplexes. Methods 2012, 57, 11−24. (65) Huppert, J. Quadruplexes in the Genome. In Quadruplex Nucleic Acids, 1st ed.; Neidle, S., Balasubramanian, S., Eds.; RSC Biomolecular Sciences; The Royal Society of Chemistry: Cambridge, United Kingdom, 2006; pp 208−227. (66) Da Silva, M. W.; Karsisiotis, A. I. Fundamentals and Applications of the Geometric Formalism of Quadruplex Folding. In Guanine Quartets: Structure and Application, 1st ed.; Fritzsche, W., Spindler, L., Eds.; The Royal Society of Chemistry: Cambridge, United Kingdom, 2013; pp 167−178. (67) Mukundan, V. T.; Phan, A. T. Bulges in G-Quadruplexes: Broadening the Definition of G-Quadruplex-Forming Sequences. J. Am. Chem. Soc. 2013, 135, 5017−5028. (68) Maizels, N. Biological Functions of G-Quadruplexes. In Guanine Quartets: Structure and Application, 1st ed.; Fritzsche, W., Spindler, L., Eds.; The Royal Society of Chemistry: Cambridge, United Kingdom, 2013; pp 215−222. (69) Biffi, G.; Tannahill, D.; McCafferty, J.; Balasubramanian, S. Quantitative visualization of DNA G-quadruplex structures in human cells. Nat. Chem. 2013, 5, 182−186. (70) Lipps, H. J.; Rhodes, D. G-quadruplex structures: in vivo evidence and function. Trends Cell Biol. 2009, 19, 414−422. (71) Russo Krauss, I.; Pica, A.; Merlino, A.; Mazzarella, L.; Sica, F. Duplex-Quadruplex Motifs in a Peculiar Structural Organization Cooperatively Contribute to Thrombin Binding of a DNA Aptamer. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2013, 69, 2403−2411. (72) Yoshida, W.; Saito, T.; Yokoyama, T.; Ferri, S.; Ikebukuro, K. Aptamer Selection Based on G4-Forming Promoter Region. PLoS One 2013, 8, e65497. 2014

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation (73) Hurley, L. H.; Wheelhouse, R. T.; Sun, D.; Kerwin, S. M.; Salazar, M.; Fedoroff, O. Y.; Han, F. X.; Han, H.; Izbicka, E.; Von Hoff, D. D. G-Quadruplexes as Targets for Drug Design. Pharmacol. Ther. 2000, 85, 141−158. (74) Rajendran, A.; Endo, M.; Hidaka, K.; Teulade-Fichou, M. P.; Mergny, J. L.; Sugiyama, H. Small molecule binding to a G-hairpin and a G-triplex: a new insight into anticancer drug design targeting G-rich regions. Chem. Commun. 2015, 51, 9181−9184. (75) Trajkovski, M.; Plavec, J. Assessing Roles of Cations in GQuadruplex-Based Nanowires by NMR. J. Phys. Chem. C 2012, 116, 23821−23825. (76) Yatsunyk, L. A.; Piétrement, O.; Albrecht, D.; Tran, P. L. T.; Renčiuk, D.; Sugiyama, H.; Arbona, J. M.; Aimé, J. P.; Mergny, J. L. Guided Assembly of Tetramolecular G-Quadruplexes. ACS Nano 2013, 7, 5701−5710. (77) Livshits, G. I.; Stern, A.; Rotem, D.; Borovok, N.; Eidelshtein, G.; Migliore, A.; Penzo, E.; Wind, S. J.; Di Felice, R.; Skourtis, S. S.; Cuevas, J. C.; Gurevich, L.; Kotlyar, A. B.; Porath, D. Long-Range Charge Transport in Single G-Quadruplex DNA Molecules. Nat. Nanotechnol. 2014, 9, 1040−1046. (78) Wang, Z. F.; Li, M. H.; Hsu, S. T. D.; Chang, T. C. Structural Basis of Sodium-Potassium Exchange of a Human Telomeric DNA Quadruplex Without Topological Conversion. Nucleic Acids Res. 2014, 42, 4723−4733. (79) Risitano, A.; Fox, K. R. Influence of Loop Size on the Stability of Intramolecular DNA Quadruplexes. Nucleic Acids Res. 2004, 32, 2598− 2606. (80) Bessi, I.; Jonker, H. R. A.; Richter, C.; Schwalbe, H. Involvement of Long-Lived Intermediate States in the Complex Folding Pathway of the Human Telomeric G-Quadruplex. Angew. Chem., Int. Ed. 2015, 54, 8444−8448. (81) Bončina, M.; Hamon, F.; Islam, B.; Teulade-Fichou, M. P.; Vesnaver, G.; Haider, S.; Lah, J. Dominant Driving Forces in Human Telomere Quadruplex Binding-Induced Structural Alterations. Biophys. J. 2015, 108, 2903−2911. (82) Gabelica, V. A pilgrim’s guide to G-quadruplex nucleic acid folding. Biochimie 2014, 105, 1−3. (83) Rajendran, A.; Endo, M.; Hidaka, K.; Sugiyama, H. Direct and Single-Molecule Visualization of the Solution-State Structures of GHairpin and G-Triplex Intermediates. Angew. Chem., Int. Ed. 2014, 53, 4107−4112. (84) You, H.; Zeng, X.; Xu, Y.; Lim, C. J.; Efremov, A. K.; Phan, A. T.; Yan, J. Dynamics and stability of polymorphic human telomeric Gquadruplex under tension. Nucleic Acids Res. 2014, 42, 8789−8795. (85) Marchand, A.; Ferreira, R.; Tateishi-Karimata, H.; Miyoshi, D.; Sugimoto, N.; Gabelica, V. Sequence and Solvent Effects on Telomeric DNA Bimolecular G-Quadruplex Folding Kinetics. J. Phys. Chem. B 2013, 117, 12391−12401. (86) 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. (87) Gomez, D.; Guédin, A.; Mergny, J. L.; Salles, B.; Riou, J. F.; Teulade-Fichou, M. P.; Calsou, P. A G-quadruplex structure within the 5′-UTR of TRF2 mRNA represses translation in human cells. Nucleic Acids Res. 2010, 38, 7187−7198. (88) Kumari, S.; Bugaut, A.; Huppert, J. L.; Balasubramanian, S. An RNA G-quadruplex in the 5 ′ UTR of the NRAS proto-oncogene modulates translation. Nat. Chem. Biol. 2007, 3, 218−221. (89) Paramasivan, S.; Rujan, I.; Bolton, P. H. Circular Dichroism of Quadruplex DNAs: Applications to Structure, Cation Effects and Ligand Binding. Methods 2007, 43, 324−331. (90) Abi-Ghanem, J.; Gabelica, V. Nucleic Acid Ion Structures in the Gas Phase. Phys. Chem. Chem. Phys. 2014, 16, 21204−21218. (91) Stadlbauer, P.; Krepl, M.; Cheatham, T. E., III; Koca, J.; Šponer, J. Structural Dynamics of Possible Late-Stage Intermediates in Folding of Quadruplex DNA Studied by Molecular Simulations. Nucleic Acids Res. 2013, 41, 7128−7143. (92) Reshetnikov, R. V.; Šponer, 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. (93) Islam, B.; Sgobba, M.; Laughton, C.; Orozco, M.; Šponer, 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. (94) Bergues-Pupo, A. E.; Arias-Gonzalez, J. R.; Morón, 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. (95) Bian, Y.; Tan, C.; Wang, J.; Sheng, Y.; Zhang, J.; Wang, W. Atomistic Picture for the Folding Pathway of a Hybrid-1 Type Human Telomeric DNA G-quadruplex. PLoS Comput. Biol. 2014, 10, e1003562. (96) Novotný, J.; Kulhánek, P.; Marek, R. Biocompatible XanthineQuadruplex Scaffold for Ion-Transporting DNA Channels. J. Phys. Chem. Lett. 2012, 3, 1788−1792. (97) Haider, S.; Neidle, S. Molecular Modeling and Simulation of GQuadruplexes and Quadruplex-Ligand Complexes. Methods Mol. Biol. 2010, 608, 17−37. (98) Haider, S.; Parkinson, G. N.; Neidle, S. Molecular dynamics and principal components analysis of human telomeric quadruplex multimers. Biophys. J. 2008, 95, 296−311. (99) Phillips, K.; Dauter, Z.; Murchie, A. I. H.; Lilley, D. M. J.; Luisi, B. J. Mol. Biol. 1997, 273, 171−182. (100) Haider, S.; Parkinson, G. N.; Neidle, S. Crystal Structure of the Potassium Form of an Oxytricha nova G-quadruplex. J. Mol. Biol. 2002, 320, 189−200. (101) 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. (102) Luu, K. N.; Phan, A. T.; Kuryavyi, V.; Lacroix, L.; Patel, D. J. Structure of the Human Telomere in K+ Solution: An Intramolecular (3 + 1) G-Quadruplex Scaffold. J. Am. Chem. Soc. 2006, 128, 9963− 9970. (103) Schultze, P.; Macaya, R. F.; Feigon, J. Three-dimensional Solution Structure of the Thrombin-binding DNA Aptamer d(GGTTGGTGTGGTTGG). J. Mol. Biol. 1994, 235, 1532−1547. (104) Hazel, P.; Parkinson, G. N.; Neidle, S. Topology Variation and Loop Structural Homology in Crystal and Simulated Structures of a Bimolecular DNA Quadruplex. J. Am. Chem. Soc. 2006, 128, 5480− 5487. (105) Schafmeister, C. E. A. F.; Ross, W. S.; Romanovski, V. LEaP; University of California: San Francisco,1995. (106) Case, D. A.; Berryman, J. T.; Betz, R. M.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Luchko, T.; Luo, R.; Madej, B.; Merz, K. M.; Monard, G.; Needham, P.; Nguyen, H.; Nguyen, H. T.; Omelyan, I.; Onufriev, A.; Roe, D. R.; Roitberg, A.; Salomon-Ferrer, R.; Simmerling, C. L.; Smith, W.; Swails, J.; Walker, R. C.; Wang, J.; Wolf, R. M.; Wu, X.; York, D. M.; Kollman, P. A. AMBER 2015; University of California: San Francisco, 2015. (107) Chung, L. W.; Hirao, H.; Li, X.; Morokuma, K. The ONIOM Method: Its Foundation and Applications to Metalloenzymes and Photobiology. WIREs Comput. Mol. Sci. 2012, 2, 327−350. (108) Vreven, T.; Byun, K. S.; Komáromi, I.; Dapprich, S.; Montgomery, J. A.; Morokuma, K.; Frisch, M. J. Combining Quantum Mechanics Methods with Molecular Mechanics Methods in ONIOM. J. Chem. Theory Comput. 2006, 2, 815−826. (109) Vreven, T.; Morokuma, K. On the Application of the IMOMO (Integrated Molecular Orbital + Molecular Orbital) Method. J. Comput. Chem. 2000, 21, 1419−1432. (110) Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; Li, H. B.; Ding, L.; Morokuma, K. The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678−5796. 2015

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016

Article

Journal of Chemical Theory and Computation (111) Tao, J. M.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 2003, 91, 146401. (112) Schaefer, A.; Huber, C.; Ahlrichs, R. Fully optimized contracted Gaussian-basis sets of triple zeta valence quality for atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−35. (113) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (114) Joung, I. S.; Cheatham, T. E., III Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020−9041. (115) 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 Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (116) Vreven, T.; Morokuma, K.; Farkas, Ö .; Schlegel, H. B.; Frisch, M. J. Geometry optimization with QM/MM, ONIOM and other combined methods. I. Microiterations and constraints. J. Comput. Chem. 2003, 24, 760−69. (117) Gkionis, K.; Platts, J. A. QM/MM studies of cisplatin complexes with DNA dimer and octamer. Comput. Theor. Chem. 2012, 993, 60−65. (118) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (119) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (120) Sinnecker, S.; Rajendran, A.; Klamt, A.; Diedenhofen, M.; Neese, F. Calculation of solvent shifts on electronic g-tensors with the conductor-like screening model (COSMO) and its self-consistent generalization to real solvents (direct COSMO-RS). J. Phys. Chem. A 2006, 110, 2235−2245. (121) Klamt, A. The COSMO and COSMO-RS solvation models. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 699−709. (122) Klamt, A.; Shüürmann, G. COSMO: A new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2 1993, 2, 799−805. (123) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic structure calculations on workstation computers: The program system Turbomole. Chem. Phys. Lett. 1989, 162, 165−169. (124) Klamt, A.; Jonas, V. Treatment of the outlying charge in continuum solvation models. J. Chem. Phys. 1996, 105, 9972−9981. (125) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal solvation model based on solute electron density and a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (126) Neese, F. ORCA-an Ab Initio, Density Functional and Semiempirical Program Package, version 2.9.1; University of Bonn: Bonn, Germany, 2012.

(127) Kruse, H.; Grimme, S. A geometrical correction for the interand intra-molecular basis set superposition error in Hartree-Fock and density functional theory calculations for large systems. J. Chem. Phys. 2012, 136, 154101. (128) Roe, D. R.; Cheatham, T. E., III PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9, 3084−3095. (129) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (130) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera-A visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605−1612. (131) Hocquet, A.; Ghomi, M. The peculiar role of cytosine in nucleoside conformational behavior: Hydrogen bond donor capacity of nucleic bases. Phys. Chem. Chem. Phys. 2000, 2, 5351−5353. (132) Hocquet, A. Intramolecular hydrogen bonding in 2′deoxyribonucleosides: an AIM topological study of the electron density. Phys. Chem. Chem. Phys. 2001, 3, 3192−3199. (133) Louit, G.; Hocquet, A.; Ghomi, M. Intramolecular hydrogen bonding in ribonucleosides: an AIM topological study of the electron density. Phys. Chem. Chem. Phys. 2002, 4, 3843−3848. (134) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999−3093. (135) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (136) Iacovelli, F.; Falconi, M. Decoding the Conformation-linked Functional Properties of Nucleic Acids by the Use of Computational Tools. FEBS J. 2015, 282, 3298−3310.

2016

DOI: 10.1021/acs.jctc.5b01025 J. Chem. Theory Comput. 2016, 12, 2000−2016