Improved Sugar Puckering Profiles for Nicotinamide Ribonucleoside

Subscriber access provided by Northern Illinois University

Article

Improved Sugar Puckering Profiles for Nicotinamide Ribonucleoside for Hybrid QM/MM Simulations Yaron Pshetitsky, Reuven Eitan, Gilit Verner, Amnon Kohen, and Dan Thomas Major J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00401 • Publication Date (Web): 04 Aug 2016 Downloaded from http://pubs.acs.org on August 5, 2016

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

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

Page 1 of 35

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

Journal of Chemical Theory and Computation

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

1

Department of Chemistry and the Lise Meitner-Minerva Center of Computational Quantum Chemistry, Bar-Ilan University, Ramat-Gan 52900, Israel 2

Department of Chemistry, University of Iowa, Iowa City, Iowa 52242 †

These authors contributed equally

[email protected] RECEIVED DATE TITLE RUNNING HEAD: Improved Sugar Puckering Profiles for Nicotinamide Ribonucleoside for Hybrid QM/MM Simulations. 1

Bar-Ilan University

2

University of Iowa

1 Environment ACS Paragon Plus


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

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 mechanicsmolecular mechanics (QM/MM) studies of NAD+/NADH dependent enzymes, the QM region is treated by semi-empirical (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.


Page 2 of 35

Page 3 of 35

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


Introduction 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 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 socalled anti- and syn-conformations,5-8 while the ψ angle is predominantly in the gauche-gauche (gg) conformation.9 The NR pseudorotation angle typically co-exists in two principle conformations: the C2’-exo-C3’-endo (or North, N) and C2’-endoC3’-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 solution.13



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

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 forcefield 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 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 ref. 18-33. To facilitate efficient, yet accurate, QM/MM MD simulations of enzymes, one may adopt a specific reaction parameter (SRP) semi-empirical (SE) approach.34 In such an approach one employs a SE Hamiltonian (e.g. AM135) with the canonical parameters replaced by a re-optimized 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


Page 4 of 35

Page 5 of 35

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


reductase (DHFR)19 and formate dehydrogenase (FDH)47 using AM1-SRP Hamiltonians. In these studies we re-optimized 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 semi-empirical 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,



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

the correction scheme is not 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 gasphase 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 conformational states defined by the above dihedral angles, further conformational sub-states, differing in their hydrogen bond (HB) patterns, were identified.


Page 6 of 35

Page 7 of 35

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


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, 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-N1C1’-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



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

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 re-optimized 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 highlevel (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 AM1-SRP levels, where the latter Hamiltonian was derived for the DHFR reaction.19 For these HL and LL methods, a two-dimensional 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


Page 8 of 35

Page 9 of 35

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


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 for the MD simulations is described by a hybrid QM(SE)/MM Hamiltonian:17

Hˆ = Hˆ QM + Hˆ MM + Hˆ QM /MM

(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 (mAM1-SRP). 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.



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

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 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 multi-dimensional 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).

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


Page 10 of 35

Page 11 of 35

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


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 gasphase results are not in nearly as good 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 M062X.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 semi-empirical 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



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

observed in our simulations, where we see a single minimum on the puckering surface in all simulations (Fig. 1 and S1), with φ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 (Fig. 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

Improved Sugar Puckering Profiles for Nicotinamide Ribonucleoside

Recommend Documents