Improved Sugar Puckering Profiles for Nicotinamide Ribonucleoside

Article pubs.acs.org/JCTC

Improved Sugar Puckering Profiles for Nicotinamide Ribonucleoside for Hybrid QM/MM Simulations Yaron Pshetitsky,†,§ Reuven Eitan,†,§ Gilit Verner,† Amnon Kohen,‡ and Dan Thomas Major*,† †

Department of Chemistry and the Lise Meitner-Minerva Center of Computational Quantum Chemistry, Bar-Ilan University, Ramat-Gan 52900, Israel ‡ Department of Chemistry, University of Iowa, Iowa City, Iowa 52242, United States S Supporting Information *

ABSTRACT: The coenzyme nicotinamide adenine dinucleotide (NAD+) and its reduced form (NADH) play ubiquitous roles as oxidizing and reducing agents in nature. The binding, and possibly the chemical redox step, of NAD+/NADH may be influenced by the cofactor conformational distribution and, in particular, by the ribose puckering of its nicotinamide− ribonucleoside (NR) moiety. In many hybrid quantum mechanics−molecular mechanics (QM/MM) studies of NAD+/NADH dependent enzymes, the QM region is treated by semiempirical (SE) methods. Recent work suggests that SE methods do not adequately describe the ring puckering in sugar molecules. In the present work we adopt an efficient and practical strategy to correct for this deficiency for NAD+/NADH. We have implemented a cost-effective correction to a SE Hamiltonian by adding a correction potential, which is defined as the difference between an accurate benchmark density functional theory (DFT) potential energy surface (PES) and the SE PES. In practice, this is implemented via a B-spline interpolation scheme for the grid-based potential energy difference surface. We find that the puckering population distributions obtained from free energy QM(SE)/MM simulations are in good agreement with DFT and in fair accord with experimental results. The corrected PES should facilitate a more accurate description of the ribose puckering in the NAD+/NADH cofactor in simulations of biological systems.

■

INTRODUCTION The coenzyme nicotinamide adenine dinucleotide (NAD+) and its reduced form (NADH), as well as their phosphorylated analogues, are of great biological importance (Scheme 1). Indeed, this redox pair plays pivotal roles in the energy processes of living cells, such as the citric acid cycle, the electron transport system, deamination of amino acids in protein metabolism, the synthesis of acetyl coenzyme A from fatty acids, and numerous other biological processes.1,2 The main catalytically relevant property of the NAD redox couple is its redox potential, and this has been studied by numerous groups.3,4 Additionally, the conformational behavior of NAD+/NADH has also been studied extensively.5−9 The conformer distribution is expected to influence ligand binding, and although not electronically coupled to the chemical step in reactions involving the cofactor, is likely to affect the donor− acceptor dynamics of the reaction via conformational vibrations. The key conformational hot spots in the © 2016 American Chemical Society

nicotinamide-ribonucleoside (NR) moiety (Scheme 2) are the glycosidic dihedral angle (χ), the C4′−C5′ exocyclic dihedral angle (ψ), and the ribose pseudorotation angle.9−11 The conformational distribution between different conformers is expected to depend both on the intrinsic conformational preference of the functional moiety (i.e., base, ribose, hydroxyl/ phosphate) and interaction between these moieties.9 It is assumed that in aqueous solution the χ angle is populated both in the so-called anti- and syn-conformations,5−8 while the ψ angle is predominantly in the gauche−gauche (gg) conformation.9 The NR pseudorotation angle typically coexists in two principle conformations: the C2′-exo-C3′-endo (or North, N) and C2′-endo-C3′-exo (or South, S) (Scheme 3).9,12 The barrier between these two conformations is in general small, and interchange occurs on a rapid time-scale in Received: April 21, 2016 Published: August 4, 2016 5179

DOI: 10.1021/acs.jctc.6b00401 J. Chem. Theory Comput. 2016, 12, 5179−5189

Article

Journal of Chemical Theory and Computation

enzymatic reactions involving the NAD redox pair, the NR moiety, along with additional parts of the active site, need to be treated by quantum mechanics (QM). This is usually implemented by electrostatically embedding the active site QM region in a molecular mechanics (MM) environment to yield a hybrid QM/MM potential energy surface (PES).17 Numerous QM/MM studies, involving the NAD redox pair, have appeared in the literature, e.g. see refs 18−33. To facilitate efficient, yet accurate, QM/MM MD simulations of enzymes, one may adopt a specific reaction parameter (SRP) semiempirical (SE) approach.34 In such an approach one employs a SE Hamiltonian (e.g., AM135) with the canonical parameters replaced by a reoptimized parameter set, which accurately describes the enzyme reaction at hand. For instance, this has been adopted for phosphoryl transfer reactions.36 We have also employed such a strategy to study numerous enzymatic reactions.19,20,37−46 In particular, we studied the reactions catalyzed by dihydrofolate reductase (DHFR)19 and formate dehydrogenase (FDH)47 using AM1-SRP Hamiltonians. In these studies we reoptimized the standard AM1 parameters to reproduce high-level ab initio and density functional theory (DFT) data for the nicotinamide moiety. In this initial work, we did not attempt to study the conformational behavior of the ribose part of the nicotinamide ribonucleoside explicitly. Recent work suggests that SE methods do not adequately describe the ring puckering in sugar molecules.48−57 This has been ascribed to inherent problems in these methods with conformational and ring-strain energies. In fact, SE methods such as AM135 and PM358,59 predict a planar ring conformation as the global puckering minimum, while the N and S conformations are not stationary points. Govender et al. presented a SRP targeting chemical glycobiology,53 and this Hamiltonian was applied to modeling monosaccharide structures, carbohydrate ring puckering, amino acid proton transfer, DNA base pair interactions, carbohydrate-aromatic πinteractions, and phosphates that are common in glycosyltransferases.54 The strength of this approach is accuracy that is similar to high-level methods, but at a fraction of the cost. However, it is challenging to derive robust semiempirical parameters for a large set of different attributes. Huang et al. adopted another approach, in their study of ribose and deoxyribose puckering in the context of RNA and DNA simulations.55 These authors introduced pucker corrections into several SE models via B-spline interpolation of the potential energy difference surface relative to benchmark ab initio QM data. The advantage of this approach is that one does not have to modify the existing Hamiltonian, which often has a fine balance of SE parameters. Hence, a localized correction is added that deals with a particular deficiency in the original method, without reducing its accuracy, and the cost of the correction is negligible. However, the correction scheme is not

Scheme 1. Reduced and Oxidized Forms of Nicotinamide Adenine Dinucleotide

Scheme 2. Reduced and Oxidized Forms of Nicotinamide− Ribonucleoside

solution.13 Interestingly, the ribose sugar puckering has subtle effect on the conformation of NR, as it causes slight changes in the relative orientation of the ribose and nicotinamide ring. Hence, sugar puckering might influence ligand binding, as well as modulate the hydride donor−acceptor distance during the reaction between the NAD+/NADH cofactor and substrate in an enzyme environment. We note that the conformational population distribution for protein-bound NAD+/NADH is likely influenced by tight interactions with the protein matrix, and might deviate significantly from that in solution. MacKerell and co-workers realized the importance of facilitating biosimulations with NAD+ and NADH and developed force field parameters for these moieties.14 Such developments have enabled extended classical molecular dynamics (MD) simulations of biological systems containing these molecules.15,16 However, to allow the treatment of

Scheme 3. Depiction of Equilibrium between North and South Conformations of the Nicotinamide−Ribonucleoside Moiety

5180


Article


and the remainder of the cyclic moiety is decomposed into triangular flaps. The puckering of the ring is defined as a combination of the angles between the flaps and the reference plane. The reference plane is determined by the coordinates of the three atoms O4′, C2′, and C3′. ϕ1 is the angle between the reference plane and the flap plane defined by O4′, C4′, and C3′, while ϕ2 is the angle between the reference plane and the flap plane defined by the atoms O4′, C1′, and C2′ (see Figure A1 and Scheme 2). This description is equivalent to using the two dihedral angles C2′−C3′−O4′−C4′ (ϕ1) and C3′−C2′− O4′−C1′ (ϕ2). Two puckering minima are expected9 when the angles ϕ1 and ϕ2 are both approximately equal to 0.5 or −0.5 radians (Scheme 3). These are noted as North, N, and South, S, conformers, respectively. We also use the notation West and East to describe the hemispheres connecting the N and the S regions of puckering space. Regarding the glycosidic dihedral angle, χ, we define this angle as C6−N1−C1′−O4′ (Scheme 2). Stable conformations are expected in the so-called syn and anti conformational regions.10 We define conformations as belonging to the anti basin if the dihedral angle is in the range 60 ± 60° and as syn if the angle is in the range 240 ± 60°. Finally, three stable conformers are expected for the exocyclic dihedral angle O4′− C4′−C5′−O5′ (ψ).10 The orientation of O5′ with respect to O4′ and C3′ can be described as gauche−gauche (gg), gauche− trans (gt), and trans−gauche (tg) (Scheme 4). NR and NRH PES Scan in the Gas-Phase and in Continuum Solvent. For each of the North and South puckering conformers the rotations around the χ and ψ dihedral angles were scanned with the Gaussian 09 ModRedundant keyword using the M06-2X functional63 in conjunction with the 6-31+G(d,p) basis set.64 Six minima were observed by the scan for each of the North and South conformers, corresponded to syn/anti states for χ and gg/gt/tg for ψ. The resulting 12 conformers for each of the NR and NRH redox states were reoptimized using the above computational method in both the gas phase and in water using the polarizable continuum model65 (PCM). The free energies were computed based on the electronic energies and standard expressions for zero-point vibrational energies and thermal corrections.66 The DFT results were also compared with MP2/ 6-311+G(2d,p) data. Reference Low and High Level Calculations for Ribose Puckering. In order to derive a PES correction surface, it is necessary to perform high-level (HL) and low-level (LL) calculations on a grid. In practice, tetrahydrofuran was employed to model the puckering of the ribose molecule. As the computational HL we adopted the M06-2X functional63 with a 6-31+G(d,p) basis set,64 to be consistent with earlier calculations on related systems.19,20 Additionally, the analogous puckering LL surfaces were computed at the AM135 and AM1SRP levels, where the latter Hamiltonian was derived for the DHFR reaction.19 For these HL and LL methods, a twodimensional puckering surface was computed with grid intervals of 5°. The difference surfaces were then computed as the difference between the DFT-based and SE-based surfaces. The model calculations used the Gaussian 09 program.67 Details regarding the implementation of the correction surface, including calculation of gradients and surface boundary treatment, may be found in the Appendix. Hybrid QM(SE)/MM Potential Energy Surface for Explicit Solvent Free Energy Simulations. The potential energy surface

particularly general, and is a patch that must be applied on a per-system basis. In the current study we adopt this latter approach due to its simplicity to present a modified SE Hamiltonian, which is applicable to protein QM/MM simulations involving NAD+/ NADH. We present free energy conformational surfaces based on multiscale QM(SE)/MM simulations in aqueous solution. The results are compared with DFT and experimental NMR data.

■

METHODOLOGY Model QM Calculations in Gas Phase and Water. Overview. In order to study the conformational behavior of NR and NRH, we adopted a combined strategy. Initially, the molecules were studied at the DFT level in the gas-phase and in aqueous solution employing a continuum solvent. Subsequently, the molecules were studied in an explicit water environment using a hybrid QM(SE)/MM strategy in conjunction with free energy MD simulations. In this latter part of the study, a puckering correction surface was added to the hybrid Hamiltonian to account for the inability of SE methods to correctly treat sugar puckering. Model Systems. We modeled the conformational behavior of NR and reduced NR (NRH) in the gas phase and in water. These molecules coexist as several different conformers in aqueous solution at room temperature.9 These conformers are described by the glycosidic dihedral angle N1−C1′ (χ), the C4′−C5′ exocyclic dihedral angle (ψ), and the ribose puckering coordinates (Schemes 2−4). In addition to the principle Scheme 4. Conformational States for the Exocyclic C4′−C5′ Dihedral Anglea

a

gg: R1 = O5′, R2 = R3 = H. tg: R2 = O5′, R1 = R3 = H. gt: R3 = O5′, R1 = R2 = H.

conformational states defined by the above dihedral angles, further conformational substates, differing in their hydrogen bond (HB) patterns, were identified. Nevertheless, since the free energy differences between these states are small, we present only the lowest free energy states of our model DFT gas-phase and continuum solvation calculations. In the hybrid QM(SE)/MM MD simulations discussed below, these are naturally included. Puckering and Definition of Conformational Coordinates. Several possible definitions of sugar puckering exist.12,55,60−62 In this work we adopt the approach of Hill and Reilly who used a triangular decomposition strategy.62 According to this approach, three selected ring atoms define a reference plane, 5181


Article


Table 1. Conformational Distribution of NRH and NR Computed at the M06-2X/6-31+G(d,p) and QM(mAM1)/MM Levels Compared to Experimental NMR Data NRH M06-2X/6-31+G(d,p) (ϕ1, ϕ2) North

South

ψ gg gg gt gt tg tg total gg gg gt gt tg tg total

NR

QM(mAM1)/MM

χ

gas-phase

water

water

anti syn anti syn anti syn

0.003 0.458 0.001 0.019 0.053 0.270 0.805 0.001 0.183 0.002 0.000 0.004 0.006 0.195

0.063 0.094 0.052 0.049 0.119 0.099 0.475 0.008 0.023 0.208 0.001 0.113 0.173 0.525

0.032 0.024 0.025 0.025 0.037 0.037 0.180 0.189 0.105 0.085 0.102 0.210 0.129 0.820

anti syn anti syn anti syn

exp

0.25

0.75

QM(mAM1)/MM

gas-phase

water

water

0.007 0.008 0.000 0.000 0.000 0.000 0.015 0.143 0.841 0.000 0.000 0.000 0.001 0.985

0.064 0.103 0.004 0.002 0.001 0.001 0.174 0.235 0.551 0.029 0.005 0.003 0.004 0.826

0.011 0.015 0.001 0.015 0.001 0.007 0.050 0.031 0.309 0.042 0.375 0.008 0.190 0.950

exp

0.49

0.51

to either of the dihedrals ϕ1 and ϕ2. Statistics were sorted into bins of width 0.01 radians for ϕ1 and ϕ2, and 0.1 radians for χ and ψ. PMF surfaces were computed using a multidimensional version of the weighted histogram analysis method (WHAM) algorithm.78,79 The free energy surfaces were depicted using the Matlab suite of programs (version 2013).

for the MD simulations is described by a hybrid QM(SE)/MM Hamiltonian:17 Ĥ = Ĥ QM + ĤMM + Ĥ QM/MM

M06-2X/6-31+G(d,p)

(1)

Four different SE methods were employed to treat the QM region containing NR and NRH: The Austin Model 135 (AM1) method, AM1-SRP, 19 AM1 with puckering correction (mAM1), and AM1-SRP with puckering correction (mAM1SRP). The MM region was modeled with the CHARMM68 force field parameters for the chlorine ion69 (for NR) and the water molecules were treated using the TIP3P model.70 Molecular Dynamics Simulations. Prior to the MD simulations, NR and NRH were solvated in a pre-equilibrated cubic water box (∼30 Å × ∼30 Å × ∼30 Å) employing periodic boundary conditions71 and TIP3P water molecules.70 The Ewald summation72 method was used to treat electrostatic interactions. In the case of oxidized NR, the system charge was neutralized by adding a chloride ion. The two systems were minimized, then gradually heated up to 298 K over the course of 25 ps, and equilibrated for 1 ns to relax the system under isothermal−isobaric (NPT) conditions using the above QM/ MM hybrid PES. The Hoover thermostat73 method was used to maintain conditions of constant pressure/temperature (CPT).74 The leapfrog integration scheme75 was employed with a time-step of 1 fs. The SHAKE algorithm76 was used to constrain all water bonds involving hydrogen atoms. Potential of Mean Force (PMF) Simulations. The free energy surfaces were obtained using umbrella sampling (US).77 In order to cover the two-dimensional puckering surface 15 US simulation windows of 200 ps per window were performed between −1 and 1 radians for both angles (ϕ1 and ϕ2). Additional windows of 1 ns each were added to thoroughly scan the low free energy paths between the Northern and Southern hemispheres. Each simulation started with a short equilibration of 2 ps before sampling. To generate the puckering free energy surfaces, a range of harmonic biasing potentials was applied to maintain sampling of ϕ1 and ϕ2 within regions of interest, while no force was applied to either of the dihedrals χ or ψ. To generate the χ−ψ free energy surfaces, harmonic biasing potentials were applied to maintain sampling of χ or ψ within regions of interest, while no force was applied

■

RESULTS AND DISCUSSION Benchmark DFT Calculations for NRH and NR. Initial benchmark calculations of the conformational preferences of NRH and NR were performed using the M06-2X functional in conjunction with the 6-31+G(d,p) basis set in the gas-phase, as well as using continuum water solvation. The results are presented in Table 1 and compared with experimental data obtained from NMR coupling constants.9 In the Supporting Information we compare the DFT results with data obtained using MP2 level of theory (Table S1 and S2). The discussion below refers to the DFT data including continuum water solvation. Experimentally, the relative population of the North conformer of NRH was determined to be 25%, whereas the South conformer was found to be 75%. In comparison, the DFT calculations obtained values of 47% and 53%, respectively, corresponding to a 0.6 kcal/mol error in the relative free energy. Additionally, experiments concluded that the gg ψ conformer of NRH constituted 47% of the total ψ population, whereas DFT predicts this to make up 19% of the populace, corresponding to a 0.8 kcal/mol error in the relative free energy. The relative population of the North conformation of NR was resolved by NMR to be 49%, while the South conformation was set at 51%. According to DFT the corresponding populations were 17% and 83%, respectively, yielding a 0.9 kcal/mol error in the relative free energy. Additionally, experiments concluded that the gg ψ conformer of NRH constituted 67% of the total ψ population, whereas DFT predicts this conformer to contribute 95%, corresponding to a 1.3 kcal/mol error in the relative free energy. Hence, DFT predicts correct conformer stability trends, with relatively small errors. We note that inclusion of continuum solvation is essential, and the gas-phase results are not in nearly as good 5182


Article

Journal of Chemical Theory and Computation agreement with experiments as when solvation effects are included. Such errors are within the range of the intrinsic errors in conformational energies expected from current density functionals, such as M06-2X.63 Importantly, these benchmark values allow us to evaluate the performance of QM(SE)/MM methods with puckering corrections, in addition to comparison with available experimental data (vide infra). QM/MM PMF Simulations of NRH and NR. Simulations employing standard semiempirical QM potentials, such as AM135 or AM1-SRP,19 are not expected to provide the correct features for the ribose pseudorotation surface. This is indeed observed in our simulations, where we see a single minimum on the puckering surface in all simulations (Figures 1 and S1), with

Figure 2. Puckering potential of mean force surface (kcal/mol) obtained for (A) NRH and (B) NR with the QM(mAM1)/MM Hamiltonian.

North and South populations of 18 and 82%, respectively, corresponding to a free energy difference of 0.9 kcal/mol. The error relative to experiment is 0.2 kcal/mol.9 The free energy barrier for the South to North transformation is 2.0 kcal/mol in the Eastern hemisphere. The QM(mAM1-SRP)/MM optimized for the DHFR reaction yields populations of 11 and 89%, corresponding to a free energy difference of 1.3 kcal/mol. The error relative to experiment is 0.6 kcal/mol.9 In this case the free energy barrier for the South to North transformation is 2.8 kcal/mol in the Western hemisphere. Hence, the relative North and South populations compare well with the experimental data of Oppenheimer and Kaplan.9 Turning to the oxidized form, NR, the agreement is not quite as good, although the general trend is correct (Figures 2B and 4). Using QM(mAM1)/MM we obtain relative North and South populations of 5 and 95%, while the QM(mAM1-SRP)/MM gives populations of 4 and 96%. These populations correspond to free energy differences of 1.8 and 2.0 kcal/mol, respectively, and the error relative to experiment is 1.7 and 1.9 kcal/mol, respectively. The respective free barriers for the South to North transformation using the two hybrid methods are 3.3 and 4.0 kcal/mol, respectively. With all methods employed, we observe a shallow local minimum in the Eastern hemisphere, which is expected to be negligible populated, in agreement with experiments for related ribose systems.80 An additional important aspect of the conformational behavior of NRH and NR is the relative orientation around the glycosidic dihedral angle (χ) and the C4′−C5′ exocyclic dihedral angle (ψ). To investigate whether the general features of these angles are correctly accounted for with our multiscale potential energy surfaces, we generated free energy surfaces using these angles as variables as well. The resulting surfaces are displayed in Figures 5 and S3−5. The only experimental

Figure 1. Puckering potential of mean force surface (kcal/mol) obtained for (A) NRH and (B) NR with the QM(AM1)/MM Hamiltonian.

ϕ1 ≈ 0.0 and ϕ2 ≈ 0.0. However, inclusion of the puckering correction potential, based on DFT calculations on tetrahydrofuran (see Appendix), in the free energy simulations rectifies this deficiency. Indeed, in all simulations a qualitatively correct puckering surface with minima in the Northern and Southern hemisphere appear (Figures 2 and S2). The North conformer is identified as a local minima with ϕ1 > 0 and ϕ2 > 0, while the South conformer is identified as a local minima with ϕ1 < 0 and ϕ2 < 0 (Figure 3 and 4). Transitions between the North and South conformations can occur via the Western or Eastern hemispheres. Examining the PMFs obtained for NRH in greater detail (Figures 2A and 3, Table 1), we may conclude that the new hybrid QM(mAM1)/MM and QM(mAM1-SRP)/MM surfaces fare as well as DFT, at a considerable lower cost per energy calculation. Employing QM(mAM1)/MM we obtain relative 5183


Article


Figure 3. Puckering in NRH from QM(mAM1)/MM MD simulations. (A) North and (B) South conformations.

information in this case is regarding the ψ angle, which shows a slight preference for the gg conformer, as mentioned above.9 Inspection of the QM(SE)/MM (SE = AM1, mAM1, AM1SRP, mAM1-SRP) NRH and NR PMF surfaces clearly show that there are distinct stationary points with interspersed local minima separated by barrier regions. These minima correspond to the syn and anti orientations around the glycosidic angle, χ, and the gg, gt, and tg states for the exocyclic dihedral angle, ψ. Specifically, the anti orientation corresponds to a χ value of ca. 0 radians, while syn corresponds to ca. ± π radians. However, these are rather wide and shallow minima, depending on the exocyclic ψ angle, as well as some dependence on the puckering, as seen by comparing the surfaces with and without puckering corrections. We note, however, that inclusion of the puckering correction does not change the qualitative features of the surfaces. Importantly, we do not observe any clear preference for either the syn or anti conformation of the glycosidic angle in all simulations of NRH and NR, consistent with experiments. Regarding the exocyclic dihedral angle, ψ, we observe minima at ca. π/3 (gg), −π/3 (gt), and π (tg), separated by slight free energy barriers. The total population of the NRH gg conformer is found to be 35% and 43% from our free energy simulations (QM(mAM1)/MM and QM(mAM1SRP)/MM, respectively), compared to 47% from experiment.9 In the case of NR the gg conformer population is 37% and 39%, respectively using the same Hamiltonians, compared to 67% from experiment.9 Seemingly, the AM1 and mAM1 based QM/

Figure 4. Puckering in NR from QM(mAM1)/MM MD simulations. (A) North and (B) South conformations.

MM Hamiltonians show a slight preference for the gt conformer, although populates all three rotamers significantly. The AM1-SRP and mAM1-SRP based QM/MM Hamiltonians show a preference for gg and gt, while somewhat less for tg. Based on the results presented, we conclude that the current approach correctly describes the main conformational features of NRH and NR. The significance of the current approach lies in providing a qualitatively correct conformational behavior for nicotinamide cofactors, whereas earlier related QM/MM studies did not consider puckering motions at all. Moreover, correct modeling of puckering motions, may facilitate a more realistic treatment of nicotinamide cofactors and their interaction with other molecules, such as enzymes and bound substrates.

■

CONCLUSIONS In this work we presented improved ribose puckering potential energy surfaces for semiempirical (SE) Hamiltonians, using DFT as a benchmark. The current results confirm earlier reports that standard SE potentials do not provide correct sugar puckering profiles. We employ a grid-based potential energy 5184


Article


Figure A1. Definition of ribose atom numbers. Dihedral ϕ1 C2′− C3′−O4′−C4′; ϕ2 C3′−C2′−O4′−C1′.

F=m

d2 r dt 2

(A1)

where m is the mass, r is the position vector of the particle, and t is the time. The force is calculated according to the gradient of the potential, V: F = −∇r V

(A2)

We can define a correction potential Ω to the LL potential surface, Vl, for a subgroup of correction parameters, X, as Ω(Χ(r)) = Vh(Χ(r)) − Vl(Χ(r))

(A3)

where Vh is the HL potential energy surface. The force in eq A2 is calculated from the composite potential, V = Vl + Ω, where Vl might be a hybrid QM(AM1)/MM potential. Hence, according to eq A2, we can define a correction force Fc to the low level method as ∂Ω ∂Ω ∂Χ =− (A4) ∂r ∂Χ ∂r Accordingly, the combined force in eq A2 is given by addition of the low level (Fl) to the correction (Fc) forces: F = Fl + Fc (A5) Fc = −

Implementing the Correction Method for MD Simulations on a Grid

Figure 5. Potential of mean force surface (kcal/mol) for the dihedral angles χ and ψ obtained for (A) NRH and (B) NR with the QM(mAM1)/MM Hamiltonian.

In order to implement the correction scheme, one needs to be able to calculate the correction force in eq A4. This requires the calculation of two quantities: the gradient of the correction potential, Ω with respect to the coordinate system, X, i.e. ( ∂Ω ),

surface correction in hybrid QM(SE)/MM free energy simulations of oxidized and reduced nicotinamide-ribonucleoside (NR/NRH) in aqueous solution. The correction surface and its gradients are calculated on the fly during the molecular dynamics runs at a minimal cost. The multiscale simulation results are in good agreement with DFT calculations with continuum water solvation and nicely accords the conformations and relative trends reported from experimental NMR data. The puckering correction introduced in this work provides qualitatively correct puckering behavior and should be useful in simulations of ligand binding and conformational effects in enzymatic reactions involving nicotinamide cofactors.

∂Χ

and the coordinate transformation from X to r, i.e. ( ∂Χ ). First, ∂r we will elaborate on how one can estimate the former quantity, followed by a discussion of the latter. Gradient of the Correction Potential

The main difficulty inherent to the correction scheme arises from the discrete nature of the numerical calculation of the grid-based correction, which is incorporated into a continuous trajectory method. During the course of a MD simulations, the system will inevitable visit spatial configuration points that are either not exactly on a grid point or possibly beyond the dimensions of the grid, i.e. X(r(τ)) ≠ X(rj), where τ represents time. Hence an interpolation scheme must be adopted. Regions beyond the dimensions of the correction potential are treated by a damping function, which smoothly switches off the correction potential. Henceforth, we will restrict the discussion to two-dimensional systems. We follow the work of Huang et al.55 and use a cubic SPLINE interpolator method. Given a function of two dimensions on a grid, f j,k  f(xj, yk), the interpolation for an off-grid point f  f(x, y) according to a two-dimensional SPLINE (2D-SPLINE) interpolation is81

■

APPENDIX In spite of the success of SRP approaches, these approximated potentials can nonetheless provide a poor description of important features of a physical system if not included in the parametrization, which in turn might lead to incorrect physical behavior and dynamics. In such a case one may improve the LL computational method, by applying an ad-hoc correction to the potential from a HL computational method as a difference potential between the LL and HL methods. In practice, such a correction scheme may be solved on a grid as a function of predefined coordinates. Such an approach for correction of an incorrect puckering PES is described below. The classical equations of motion is governed by the force, F, according to Newton’s second law:

f = Ax f j + Bx f j + 1 + Cxf jxx + Dx f jxx+ 1

(A6)

where 5185


Article

Journal of Chemical Theory and Computation f j = A y f j , k + By f j + 1, k + Cyf jyy, k + Dyf jyy+ 1, k

Coordinate Transformation

(A7)

The coordinate transformation from X to r ( ∂Χ ) in eq A4 is ∂r more complicated and depends on the specific coordinates employed (i.e. X), and of the system at hand. Here, we will present two relevant parameters, X, which exemplify how this grid correction scheme may be used. In the following, X is defined as the puckering coordinates of a five-membered ring. Here, we use the Hill and Reilly approach,62 where the puckering is defined by two dihedral angles. The coordinate system, X, is defined by two angles: ϕ1, the dihedral angle between atoms C2′−C3′−O4′−C4′, and ϕ2, the dihedral angle between atoms C3′−C2′−O4′−C1′ (see Figure A1). The dihedral angle between ri, rj, rk, and rl (ϕijkl) is formally defined as82

and Ax =

xj + 1 − x xj + 1 − xj

Bx = 1 − Ax ;

;

1 (Ax 3 − Ax )(xj + 1 − xj)2 ; 6 1 Dx = (Bx 3 − Bx )(xj + 1 − xj)2 6 yj + 1 − y ; By = 1 − A y ; Ay = yj + 1 − yj Cx =

1 (A y 3 − A y )(yj + 1 − yj )2 ; 6 1 Dy = (By 3 − By )(yj + 1 − yj )2 6

(A8a)

Cy =

ϕijkl = cos−1(t ·̂ û )

where t ̂ and û are the unit vectors of t = rij × rkj

(A8b)

yy f xx j and f k,j are the numerical second order derivatives with respect to the relevant interpolated axis. These elements are calculated numerically by a SPLINE interpolation subroutine. The derivative of f with respect to x is straightforward:

3Bx 2 − 1 (xj + 1 − xj)f jxx+ 1 6

where ∂f j ∂y

∂f j ∂y

=

(A9)

yk + 1 − yk −

∂f jxx ∂y

∂r

Avoiding Discontinuities in the Forces at the Correction Potential Boundaries

(A10)

The correction scheme is only applied to certain regions of the configuration space of X. For a continuous potential function, the correction function should decay smoothly towards zero at the edges of the correction subspace. Additionally, one should ensure that the damping function also maintains continuity of the forces. We choose a damping function of the form:

is fj,k+1 − fj,k 3By2 − 1 6

−

3A y2 − 1 6

(yk + 1 − yk )f jyy, k

(yk + 1 − yk )f jyy, k + 1

⎧ 1, r∈S ⎪ D(r) = ⎨ ⎪ −(κ(r))2 , r∉S ⎩e

(A11)

has to be evaluated numerically. Nevertheless, we would

Λ(r) = D(r)Ω(r)

∂f j ∂y

(A17)

and we get the desired continuity requirements; Ω and D, and hence Λ, are continuous funnctions. The derivative of Λ is ⎧ ∂Ω , r∈S ⎪ ⎪ ∂r ∂Λ =⎨ ⎞ ∂r ∂κ ⎪ −(κ(r))2⎜⎛ ∂Ω − 2κ(r) Ω(r)⎟ , r ∉ S ⎪e ⎝ ⎠ ∂r ∂r ⎩

∂f jxx

The function

(A16)

where S contains all the points inside a close region σ(r), and κ(σ(r)) = 0. We now redefined the correction potential as

like to minimize the numerical contribution to this derivative in order to avoid unnecessary numerical errors. Therefore, we do not calculate the numerical first derivative of f xx j with respect to y directly. As explained above, f xx is a numerical second j derivative given by the SPLINE interpolation scheme. Hence, deriving this quantity numerically will result in a third order numerical error. Instead, we use the fact that x and y are independent, and we can exchange the derivation order: ⎛ ∂f j ⎞xx = ⎜⎜ ⎟⎟ ∂y ⎝ ∂y ⎠

(A15)

where rnm = rn − rm. On initial inspection, one might presume that the cos−1 function, the derivative of ϕijkl with respect to r, should reach infinity at ϕijkl = kπ (k is an integer number). This is not the caseclearly the ring can flip and the forces should not go to infinity. Swope and Ferguson have derived a formalism that overcomes this difficulty and facilitates the calculation of these derivatives.82 We therefore have used eq 40−43 of ref 82 to calculate the terms of ∂Χ .

The derivative with respect to y is given by ∂f j ∂f j + 1 ∂f jxx ∂f jxx+ 1 ∂f = Ax + Bx + Cx + Dx ∂y ∂y ∂y ∂y ∂y

(A14)

u = rlk × rkj

fj+1 − fj 3Ax 2 − 1 ∂f = − (xj + 1 − xj)f jxx ∂x xj + 1 − xj 6 −

(A13)

(A12)

Hence,

is given explicitly in eq A11. We can therefore

lim

perform the second derive numerically using the second numerical derivative given by the SPLINE interpolation. Now, the numerical error will be of second order rather than third order, as would be the case if we directly perform the derivation of f xx j numerically with respect to y.

∂Λ ∂r

r → σ(r)−

(A18)

is also continuous, where at r = σ(r): ∂Λ ∂Ω ∂Λ = = lim + r → σ(r) ∂r ∂r ∂r

(A19) −

where we use the fact that κ(σ(r)) = 0 (σ(r) is inside the shell, i.e. σ(r)− ∈ S, and σ(r)+ is outside the shell, i.e. σ(r)+ ∉ S). In practice we have implemented the damping function as 5186


Article

Journal of Chemical Theory and Computation ⎧ 1, r ≤ r0 ⎪ D(r) = ⎨ 2 ⎪ −(r − r0) , r > r0 ⎩e

Insights into the Mechanism of Bovine CD38/NAD+Glycohydrolase from the X-Ray Structures of Its Michaelis Complex and CovalentlyTrapped Intermediates. PLoS One 2012, 7, e34918. (11) Saenger, W. Principles of Nucleic Acid Structure; Springer-Verlag: New York, 1984. (12) Altona, C.; Sundaralingam, M. Conformational analysis of the sugar ring in nucleosides and nucleotides. New description using the concept of pseudorotation. J. Am. Chem. Soc. 1972, 94, 8205−8212. (13) Altona, C.; Sundaralingam, M. Conformational analysis of the sugar ring in nucleosides and nucleotides. Improved method for the interpretation of proton magnetic resonance coupling constants. J. Am. Chem. Soc. 1973, 95, 2333−2344. (14) Pavelites, J. J.; Gao, J.; Bash, P. A.; Mackerell, A. D. A Molecular Mechanics Force Field for NAD+, NADH, and the Pyrophosphate Groups of Nucleotides. J. Comput. Chem. 1997, 18, 221−239. (15) Rod, T. H.; Radkiewicz, J. L.; Brooks, C. L., III Correlated motion and the effect of distal mutations in dihydrofolate reductase. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 6980−6985. (16) Khavrutskii, I.; Price, D.; Lee, J.; Brooks, C. L., III Conformational change of the methionine 20 loop of Escherichia coli dihydrofolate reductase modulates pKa of the bound dihydrofolate. Protein Sci. 2007, 16, 1087−1100. (17) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227−249. (18) Garcia-Viloca, M.; Truhlar, D. G.; Gao, J. Reaction-Path Energetics and Kinetics of the Hydride Transfer Reaction Catalyzed by Dihydrofolate Reductase. Biochemistry 2003, 42, 13558−13575. (19) Doron, D.; Major, D. T.; Kohen, A.; Thiel, W.; Wu, X. Hybrid Quantum and Classical Simulations of the Dihydrofolate Reductase Catalyzed Hydride Transfer Reaction on an Accurate Semi-Empirical Potential Energy Surface. J. Chem. Theory Comput. 2011, 7, 3420− 3437. (20) Vardi-Kilshtain, A.; Major, D. T.; Kohen, A.; Engel, H.; Doron, D. Hybrid Quantum and Classical Simulations of the Formate Dehydrogenase Catalyzed Hydride Transfer Reaction on an Accurate Semiempirical Potential Energy Surface. J. Chem. Theory Comput. 2012, 8, 4786−4796. (21) Castillo, R.; Andres, J.; Moliner, V. Catalytic Mechanism of Dihydrofolate Reductase Enzyme. A Combined Quantum-MechanicalMolecular-Mechanical Characterization of Transition State Structure for the Hydride Transfer Step. J. Am. Chem. Soc. 1999, 121, 12140− 12147. (22) Thorpe, I. F.; Brooks, C. L. Barriers to Hydride Transfer in Wild Type and Mutant Dihydrofolate Reductase from E. coli. J. Phys. Chem. B 2003, 107, 14042−14051. (23) Liu, H.; Warshel, A. The Catalytic Effect of Dihydrofolate Reductase and Its Mutants Is Determined by Reorganization Energies. Biochemistry 2007, 46, 6011−6025. (24) Luk, L. Y. P.; Ruiz-Pernía, J. J.; Dawson, W. M.; Roca, M.; Loveridge, E. J.; Glowacki, D. R.; Harvey, J. N.; Mulholland, A. J.; Tuñoń , I.; Moliner, V.; Allemann, R. K. Unraveling the role of protein dynamics in dihydrofolate reductase catalysis. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 16344−16349. (25) Ruiz-Pernia, J. J.; Luk, L. Y. P.; García-Meseguer, R.; Martí, S.; Loveridge, E. J.; Tuñoń , I. a.; Moliner, V.; Allemann, R. K. Increased Dynamic Effects in a Catalytically Compromised Variant of Escherichia coli Dihydrofolate Reductase. J. Am. Chem. Soc. 2013, 135, 18689−18696. (26) Luk, L. Y. P.; Ruiz-Pernía, J. J.; Dawson, W. M.; Loveridge, E. J.; Tuñoń , I. a.; Moliner, V.; Allemann, R. K. Protein Isotope Effects in Dihydrofolate Reductase From Geobacillus stearothermophilus Show Entropic−Enthalpic Compensatory Effects on the Rate Constant. J. Am. Chem. Soc. 2014, 136, 17317−17323. (27) Luk, L. Y. P.; Ruiz-Pernía, J. J.; Adesina, A. S.; Loveridge, E. J.; Tuñoń , I.; Moliner, V.; Allemann, R. K. Chemical ligation and isotope labeling to locate dynamic effects during catalysis by dihydrofolate reductase. Angew. Chem., Int. Ed. 2015, 54, 9016−9020.

(A20)

where r = ∥r∥, with a range defined by some radius r0.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00401. Puckering potential of mean surfaces for NRH and NR using QM(AM1-SRP)/MM and QM(mAM1-SRP)/MM Hamiltonians, potential of mean force surfaces for dihedral angles χ and ψ. QM(AM1)/MM, QM(AM1SRP)/MM, and QM(mAM1-SRP)/MM Hamiltonians, comparison of DFT and MP2 calculations, and populations computed with the QM(mAM1-SRP)/MM Hamiltonian (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions §

Y.P. and R.E. contributed equally.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work has been supported by the Israel Science Foundation (Grant no. 2146/15) and the United States−Israel Binational Science Foundation (Grant no. 2012340).

■

REFERENCES

(1) Alberts, B.; Bray, D.; Lewis, J.; Raff, M.; Roberts, K.; Watson, J. D. Molecular Biology of the Cell, 2nd ed.; Garland Publishing: New York, 1989. (2) Kroschwitz, J. L.; Winokur, M. Chemistry: General, Organic, Biological; McGraw-Hill: New York, 1985. (3) Gebicki, J.; Marcinek, A.; Zielonka, J. Transient Species in the Stepwise Interconversion of NADH and NAD+. Acc. Chem. Res. 2004, 37, 379−386. (4) Frey, P. A.; Hegeman, A. D. Enzymatic Reaction Mechanisms; Oxford University Press: New York, 2007. (5) Zens, A. P.; Williams, T. J.; Wisowaty, J. C.; Fisher, R. R.; Dunlap, R. B.; Bryson, T. A.; Ellis, P. D. Nuclear Magnetic-Resonance Studies on Pyridine Dinucleotides 0.2. Solution Conformational Dynamics of Nicotinamide Adenine-Dinucleotide and Nicotinamide Mononucleotide as Viewed by Proton T1Measurements. J. Am. Chem. Soc. 1975, 97, 2850−2857. (6) Egan, W.; Forsen, S.; Jacobus, J. The solution conformation of nicotinamide mononucleotide: a quantitative application of the nuclear Overhauser effect. Biochemistry 1975, 14, 735−742. (7) Birdsall, B.; Birdsall, N. J.; Feeney, J.; Thornton, J. A nuclear magnetic resonance investigation of the conformation of nicotinamide mononucleotide in aqueous solution. J. Am. Chem. Soc. 1975, 97, 2845−2850. (8) Lee, C.-Y.; Raszka, M. J. Determination of solution structure of diphosphopyridine coenzymes with paramagnetic shift and broadening reagents. J. Magn. Reson. 1975, 17, 151−160. (9) Oppenheimer, N. J.; Kaplan, N. O. Proton magnetic resonance study of the intramolecular association and conformation of the alpha and beta pyridine mononucleotides and nucleosides. Biochemistry 1976, 15, 3981−3989. (10) Egea, P. F.; Muller-Steffner, H.; Kuhn, I.; Cakir-Kiefer, C.; Oppenheimer, N. J.; Stroud, R. M.; Kellenberger, E.; Schuber, F. 5187


Article

Journal of Chemical Theory and Computation (28) Wang, Z.; Antoniou, D.; Schwartz, S. D.; Schramm, V. L. Hydride Transfer In Dhfr By Transition Path Sampling, Kinetic Isotope Effects, And Heavy Enzyme Studies. Biochemistry 2016, 55, 157−166. (29) Dametto, M.; Antoniou, D.; Schwartz, S. D. Barrier crossing in dihydrofolate reductase does not involve a rate-promoting vibration. Mol. Phys. 2012, 110, 531−536. (30) Agarwal, P. K.; Billeter, S. R.; Hammes-Schiffer, S. Nuclear quantum effects and enzyme dynamics in dihydrofolate reductase catalysis. J. Phys. Chem. B 2002, 106, 3283−3293. (31) Agarwal, P. K.; Billeter, S. R.; Rajagopalan, P. T.; Benkovic, S. J.; Hammes-Schiffer, S. Network of coupled promoting motions in enzyme catalysis. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 2794−2799. (32) Watney, J. B.; Agarwal, P. K.; Hammes-Schiffer, S. Effect of Mutation on Enzyme Motion in Dihydrofolate Reductase. J. Am. Chem. Soc. 2003, 125, 3745−3750. (33) Liu, C. T.; Francis, K.; Layfield, J. P.; Huang, X.; HammesSchiffer, S.; Kohen, A.; Benkovic, S. J. Escherichia coli dihydrofolate reductase catalyzed proton and hydride transfers: Temporal order and the roles of Asp27 and Tyr100. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 18231−18236. (34) Rossi, I.; Truhlar, D. G. Parameterization of NDDO wavefunctions using genetic algorithms. An evolutionary approach to parameterizing potential energy surfaces and direct dynamics calculations for organic reactions. Chem. Phys. Lett. 1995, 233, 231− 236. (35) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. AM1: A New General Purpose Quantum Mechanical Molecular Model. J. Am. Chem. Soc. 1985, 107, 3902−3909. (36) Nam, K.; Cui, Q.; Gao, J.; York, D. M. Specific Reaction Parametrization of the AM1/d Hamiltonian for Phosphoryl Transfer Reactions: H, O, and P Atoms. J. Chem. Theory Comput. 2007, 3, 486− 504. (37) Major, D. T.; York, D. M.; Gao, J. L. Solvent polarization and kinetic isotope effects in nitroethane deprotonation and implications to the nitroalkane oxidase reaction. J. Am. Chem. Soc. 2005, 127, 16374−16375. (38) Major, D. T.; Gao, J. A combined quantum mechanical and molecular mechanical study of the reaction mechanism and alphaamino acidity in alanine racemase. J. Am. Chem. Soc. 2006, 128, 16345−16357. (39) Major, D. T.; Nam, K.; Gao, J. Transition State Stabilization and [alpha]-Amino Carbon Acidity in Alanine Racemase. J. Am. Chem. Soc. 2006, 128, 8114−8115. (40) Major, D. T.; Heroux, A.; Orville, A. M.; Valley, M. P.; Fitzpatrick, P. F.; Gao, J. Differential quantum tunneling contributions in nitroalkane oxidase catalyzed and the uncatalyzed proton transfer reaction. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 20734−20739. (41) Rubinstein, A.; Major, D. T. Catalyzing Racemizations in the Absence of a Cofactor: The Reaction Mechanism in Proline Racemase. J. Am. Chem. Soc. 2009, 131, 8513−8521. (42) Doron, D.; Weitman, M.; Vardi-Kilshtain, A.; Azuri, A.; Engel, H.; Major, D. T. Multiscale Quantum-Classical Simulations of Enzymes. Isr. J. Chem. 2014, 54, 1108−1117. (43) Carvalho, A. T. P.; Barrozo, A.; Doron, D.; Vardi Kilshtain, A.; Major, D. T.; Kamerlin, S. C. L. Challenges in computational studies of enzyme structure, function and dynamics. J. Mol. Graphics Modell. 2014, 54, 62−79. (44) Engel, H.; Eitan, R.; Azuri, A.; Major, D. T. Nuclear Quantum Effects in Chemical Reactions via Higher-Order Path-Integral Calculations. Chem. Phys. 2015, 450−451, 95−101. (45) Nitoker, N.; Major, D. T. Understanding the Reaction Mechanism and Intermediate Stabilization in Serine Racemase Using Multiscale Quantum-Classical Simulations. Biochemistry 2015, 54, 516−527. (46) Vardi-Kilshtain, A.; Nitoker, N.; Major, D. T. Nuclear Quantum Effects and Kinetic Isotope Effects in Enzyme Reactions. Arch. Biochem. Biophys. 2015, 582, 18−27.

(47) Vardi-Kilshtain, A.; Azuri, A.; Major, D. T. Path-Integral Calculations of Heavy Atom Kinetic Isotope Effects in Condensed Phase Reactions Using Gradient-Based Forward Corrector Algorithms. J. Comput. Chem. 2012, 33, 435−441. (48) McNamara, J. P.; Muslim, A.-M.; Abdel-Aal, H.; Wang, H.; Mohr, M.; Hillier, I. H.; Bryce, R. A. Towards a quantum mechanical force field for carbohydrates: a reparametrized semi-empirical MO approach. Chem. Phys. Lett. 2004, 394, 429−436. (49) Barnett, C. B.; Naidoo, K. J. Stereoelectronic and Solvation Effects Determine Hydroxymethyl Conformational Preferences in Monosaccharides. J. Phys. Chem. B 2008, 112, 15450−15459. (50) Sattelle, B. M.; Almond, A. Less is more when simulating unsulfated glycosaminoglycan 3D-structure: Comparison of GLYCAM06/TIP3P, PM3CARB1/TIP3P, and SCC-DFTB-D/TIP3P predictions with experiment. J. Comput. Chem. 2010, 31, 2932−2947. (51) Barnett, C. B.; Naidoo, K. J. Ring puckering: a metric for evaluating the accuracy of AM1, PM3, PM3CARB-1, and SCC-DFTB carbohydrate QM/MM simulations. J. Phys. Chem. B 2010, 114, 17142−17154. (52) Islam, S. M.; Roy, P.-N. Performance of the SCC-DFTB Model for Description of Five-Membered Ring Carbohydrate Conformations: Comparison to Force Fields, High-Level Electronic Structure Methods, and Experiment. J. Chem. Theory Comput. 2012, 8, 2412− 2423. (53) Govender, K.; Gao, J.; Naidoo, K. J. AM1/d-CB1: A Semiempirical Model for QM/MM Simulations of Chemical Glycobiology Systems. J. Chem. Theory Comput. 2014, 10, 4694−4707. (54) Govender, K. K.; Naidoo, K. J. Evaluating AM1/d-CB1 for Chemical Glycobiology QM/MM Simulations. J. Chem. Theory Comput. 2014, 10, 4708−4717. (55) Huang, M.; Giese, T. J.; Lee, T. S.; York, D. M. Improvement of DNA and RNA Sugar Pucker Profiles from Semiempirical Quantum Methods. J. Chem. Theory Comput. 2014, 10, 1538−1545. (56) Mayes, H. B.; Broadbelt, L. J.; Beckham, G. T. How sugars pucker: electronic structure calculations map the kinetic landscape of five biologically paramount monosaccharides and their implications for enzymatic catalysis. J. Am. Chem. Soc. 2014, 136, 1008−1022. (57) Huang, M.; Giese, T. J.; York, D. M. Nucleic Acid Reactivity: Challenges for Next-Generation Semiempirical Quantum Models. J. Comput. Chem. 2015, 36, 1370−1389. (58) Stewart, J. J. P. Optimization of parameters for semiempirical methods I. Method. J. Comput. Chem. 1989, 10, 209−220. (59) Stewart, J. J. P. Optimization of parameters for semiempirical methods II. Applications. J. Comput. Chem. 1989, 10, 221−264. (60) Cremer, D.; Pople, J. A. General definition of ring puckering coordinates. J. Am. Chem. Soc. 1975, 97, 1354−1358. (61) Sato, T. Another Method for Specifying Furanose Ring Puckering. Nucleic Acids Res. 1983, 11, 4933−4938. (62) Hill, A. D.; Reilly, P. J. Puckering coordinates of monocyclic rings by triangular decomposition. J. Chem. Inf. Model. 2007, 47, 1031−1035. (63) 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. (64) Hehre, W. J. Ab initio molecular orbital theory. Acc. Chem. Res. 1976, 9, 399−406. (65) Miertuš, S.; Scrocco, E.; Tomasi, J. Electrostatic interaction of a solute with a continuum. A direct utilizaion of AB initio molecular potentials for the prevision of solvent effects. Chem. Phys. 1981, 55, 117−129. (66) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, CA, 2000. (67) Frisch, M. J.; T, 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, 5188


Article

Journal of Chemical Theory and Computation 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, J. M.; 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, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision A.2; Gaussian, Inc.: Wallingford CT, 2009. (68) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102, 3586−616. (69) Roux, B. The molecular basis of the valence selectivity of the gramicidin channel: A molecular dynamics free energy perturbation study. Biophys. J. 1996, 71, 3177−3185. (70) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (71) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: 1989. (72) Nam, K.; Gao, J.; York, D. M. An efficient linear-scaling Ewald method for long-range electrostatic interactions in combined QM/ MM calculations. J. Chem. Theory Comput. 2005, 1, 2−13. (73) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (74) Andersen, H. C. Molecular-Dynamics Simulations at Constant Pressure and-or Temperature. J. Chem. Phys. 1980, 72, 2384−2393. (75) Hockney, R. W. The Potential Calculation and Some Applications. Methods Comput. Phys. 1970, 9, 136−210. (76) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327−341. (77) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187−199. (78) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The weighted histogram analysis method for freeenergy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13, 1011−1021. (79) Doron, D.; Kohen, A.; Major, D. T. Collective Reaction Coordinate for Hybrid Quantum and Molecular Mechanics Simulations: A Case Study of the Hydride Transfer in Dihydrofolate Reductase. J. Chem. Theory Comput. 2012, 8, 2484−2496. (80) Olson, W. K. Three-state models of furanose pseudorotation. Nucleic Acids Res. 1981, 9, 1251−1262. (81) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing (1986−1992), Third ed.; Cambridg University Press: New York, NY, 2007. (82) Swope, W. C.; Ferguson, D. M. Alternative Expressions for Energies and Forces Due to Angle Bending and Torsional Energy. J. Comput. Chem. 1992, 13, 585−594.

5189


Improved Sugar Puckering Profiles for Nicotinamide Ribonucleoside

Recommend Documents