MM Approaches To Elucidate ... - ACS Publications

Article pubs.acs.org/JCTC

Different QM/MM Approaches To Elucidate Enzymatic Reactions: Case Study on ppGalNAcT2 Pavel Janoš,†,‡ Tomás ̌ Trnka,†,‡ Stanislav Kozmon,†,§ Igor Tvaroška,†,§ and Jaroslav Koča*,†,‡ †

Central European Institute of Technology (CEITEC) and ‡Faculty of Science-National Centre for Biomolecular Research, Masaryk University, 62500 Brno, Czech Republic § Institute of Chemistry, Slovak Academy of Sciences, 845 38 Bratislava, Slovak Republic S Supporting Information *

ABSTRACT: Hybrid QM/MM computational studies can provide invaluable insight into the mechanisms of enzymatic reactions that can be exploited for rational drug design. Various approaches are available for such studies. However, their strengths and weaknesses may not be immediately apparent. Using the retaining glycosyltransferase ppGalNAcT2 as a case study, we compare different methodologies used to obtain reaction paths and transition state information. Ab Initio MD using CPMD coupled with the String Method is used to derive the minimum free energy reaction path. The geometrical features of the free energy path, especially around the transition state, agree with the minimum potential energy path obtained by the much less computationally expensive Nudged Elastic Band method. The barrier energy, however, differs by 8 kcal/mol. The free energy surface generated by metadynamics provides a rough overview of the reaction and can confirm the physical relevance of optimized paths or provide an initial guess for path optimization methods. Calculations of enzymatic reactions are usually performed at best at the DFT level of theory. A comparison of widely used functionals with high-level DLPNO-CCSD(T)/CBS data on the potential energy profile serves as a validation of the usability of DFT for this type of enzymatic reaction. The M06-2X meta-hybrid functional in particular matches the DLPNO-CCSD(T)/CBS reference extremely well with errors within 1 kcal/mol. However, even pure-GGA functional OPBE provides sufficient accuracy for this type of study.

■

INTRODUCTION

A number of computational studies have been performed on glycosyltransferases to elucidate the mechanism of the reaction.7 One option is to explore the potential energy surface (PES) using static structures. This approach requires a prior definition of one to three geometrical parameters (scan coordinates) that serve as the independent variables defining the location of the system on the PES. All remaining degrees of freedom are relaxed by geometry optimization at each scan point. Although this method is conceptually simple, its success or failure is highly dependent on a suitable choice of scan coordinates. If these coordinates do not describe the reactive motions adequately, one commonly obtains a flawed (discontinuous) potential energy profile. Unfortunately, the presence of such discontinuities is often very hard to detect. One may alternatively focus only on locating the minimum potential energy reaction path using, for example, the Nudged Elastic Band method (NEB)8 or string optimization methods.9 This approach is not limited by the requirement to choose only a few geometrical parameters to describe the reaction and can be carried out easily in full Cartesian coordinate space. This ensures that a continuous reaction path will be obtained. However, there can be multiple alternative pathways linking the

Glycosyltransferases are a vast group of enzymes that participate in the process of glycan biosynthesis by catalyzing the formation of the glycosidic bond. They facilitate the transfer of a sugar moiety from activated donor substrates (sugarnucleotide) to an acceptor. The disruption of the glycan composition of cells has been linked to various diseases. Changes in O-GlcNAcylation facilitated by O-GlcNAc transferase (OGT), for example, have been associated with cancer, diabetic complications, and even neurodegenerative diseases including Alzheimer’s disease.1 O-GalNAcylation has been shown to play a role in cancer metastasis.2 The increase in branched N-glycans has been implicated in poor cancer prognosis.3 Glycosyltransferases involved in the biosynthesis of the Mycobacterium cell wall are interesting targets for antituberculosis drug development.4,5 Glycosylation is also the mechanism behind bacterial toxins produced by the Clostridium and Legionella families, the latter being the causative agent of Legionnaires’ disease.6 This places glycosyltransferases into the spotlight of drug development. Detailed insight into the reaction mechanism is needed for rational drug design. An ideal starting point for rational drug design is the structure of the transition state of the reaction. Computational methods can provide such information. © 2016 American Chemical Society

Received: May 23, 2016 Published: October 27, 2016 6062

DOI: 10.1021/acs.jctc.6b00531 J. Chem. Theory Comput. 2016, 12, 6062−6076

Article

Journal of Chemical Theory and Computation

Unfortunately, the improvement in the reliability provided by running AIMD simulations comes at an enormous computational cost compared to PES-based approaches. As a rule of thumb, one may expect that an FES-based study of a given enzymatic system will require 1 to 2 orders of magnitude more CPU time than a similar PES-based study. In this work, we have set out to determine whether this additional cost is justified by additional insight into the mechanism of the reaction. Regardless of the particular method used to study the reaction mechanism, the level of theory is usually restricted to GGA DFT. AIMD simulations, in particular, cannot be realistically performed on enzymatic QM/MM systems with hybrid density functionals or correlated methods. It is thus important to know whether the resulting reaction mechanism is affected by possible inaccuracies or artifacts of the approximate density functional used. Fortunately, the recent introduction of near-linear-scaling coupled clusters methods18,19 has enabled CCSD(T) calculations on systems containing several hundreds of atoms. This enables us to obtain a high-level potential energy profile by single-point calculations on the optimized reaction path. We will subsequently validate the more approximate DFT results against this reference profile. In this study, we have investigated the performance of several QM/MM methods in predicting the catalytic mechanism of glycosyltransferases. We have selected the polypeptide UDPGalNAc transferase (ppGalNAcT2) in its human isoform 2 as the model system for our study. The glycosyltransferase ppGalNAcT2 catalyzes the transfer of an N-acetylgalactosaminyl (GalNAc) residue from UDPGalNAc to the serine or threonine hydroxyl moieties of an acceptor protein. ppGalNAcT2 is a retaining metal-dependent enzyme and has already been subject to extensive experimental20−27 and computational investigation28−30 by several research groups. According to these results, this enzyme catalyzes a same-face nucleophilic substitution, and the reaction mechanism is very similar to other configuration-retaining glycosyltransferases. It is thus a suitable representative of this important group of enzymes for our methodological study.

reactant and product basins of the potential energy surface (PES), and there is absolutely no guarantee that the optimized reaction path is the most representative or physically significant one. To assess the validity of the optimized path, one usually needs some information about the overall shape of the PES. Thus, PES scan and path optimization methods need to be combined to ensure meaningful results. There is, however, a significant shortcoming shared by all of these methods. Because they are based on “static” structures (generated purely by a local optimization of geometry) without any conformational sampling, they are in danger of being locked into an incorrect local minimum and therefore misrepresenting the enzymatic reaction. Additionally, focusing only on the potential energy instead of the free energy neglects any entropic contributions that might alter both the observed reaction path and its barrier height. To switch to exploring the free energy surface, including the entropic contribution, Ab Initio MD (AIMD) simulations have to be performed. Due to their enormous computational cost, these are usually carried out at a semiempirical level of theory. Semiempirical methods are significantly faster than ab initio QM methods. However, the reliability of semiempirical methods in calculations of carbohydrate conformations is insufficient, though recently the PM3-CARB1 semiempirical method was developed for carbohydrates.10 While it is an improvement over previous semiempirical methods, its performance is still lacking.11,12 More recently another carbohydrate specific semiempirical method was developed, AM1/d-CB1.13 The quality of this parametrization, however, remains to be proven. As with the PES, there are various approaches for exploring the free energy surface (FES). Improved sampling methods are routinely used to overcome the time scale limitations of ab initio MD and obtain an overview of the free energy landscape. A particularly common choice in this regard is the metadynamics14 method. However, all these FES exploration approaches share the pitfalls of PES scans in terms of the suitability of the collective variables selected to describe the reaction. Minimum free energy reaction paths (MFEP) can be optimized using similar methods to their potential energy counterparts, only applied to a set of collective variables instead of full Cartesian coordinates. However, in striking contrast to FES exploration methods such as metadynamics, MFEP optimization methods do not place any limitation on the number of collective variables employed. Several tens of collective variables may thus be easily used at the same time to make their specific selection not crucial. We have previously applied the string method (STM),9,15 an MFEP optimization method, to study the reaction mechanism of OGT glycosyltransferase.16 The last group of commonly used approaches is represented by the Transition Path Sampling (TPS) method.17 Similarly to MFEP optimization, this approach focuses on sampling a small subspace of the whole free energy landscape accurately instead of trying to obtain an overview. However, instead of trying to find the path connecting reactants and products, TPS optimizes the dividing surface between the reactant and product basins. Both MFEP and TPS thus ultimately produce the same information about the approximate location of the transition state but differ in what kind of additional insight into the reaction they offer (path through the basins vs shape of the barrier).

■

METHODOLOGY System Preparation. The initial structural model for the CPMD MTD/STM taken from our previous work30 was constructed from the X-ray structure of ppGalNAcT human isoform 220 (PDB: 2FFU) containing acceptor peptide EA2 and the UDP part of the donor molecule and the X-ray structure of human isoform 1022 containing UDP-GalNAc (PDB: 2D7I). The resulting structure was protonated and contained both acceptor peptide EA2 and donor substrate UDP-GalNAc. The system was solvated in a cubic water box so that the distance between the protein and box wall was 10 Å. Nine Cl− ions were added to ensure the neutrality of the simulation box. The system was described using the FF99SB force field31 and the GLYCAM06 force field32 parameters for the GalNAc ring. Water was described using the TIP3P water model.33 The system was equilibrated following our standard protocol. Solvent molecule optimization was done using 3000 steps of steepest descent minimization; the rest of the system was held with a 50 kcal/mol restraint. The simulation box was heated to 300 K during the 100 ps-long constant volume simulation, followed by a 300 ps-long 1-bar constant pressure simulation, then a series of 10 ps-long constant pressure simulations with gradually decreasing strength of the restraint, from 25 to 10 kcal/mol, and the simulations with 5 and 2.5 6063


Article


switching to the QM/MM description. After that, the system was heated to 300 K using a Berendsen weak coupling thermostat over the course of 1.2 ps. An additional 3 ps of NVT simulation were done to equilibrate the QM/MM system using three Nose-Hoover thermostats, each coupled to one of the regions: QM zone, protein, solvent. The same thermostat setup was used in the production run. The QM/MM system equilibrated in this way served as a starting point for both the MTD and STM simulations using the same conditions. Collective Variables. During the ppGalNAcT2 enzymatic reaction, the phosphoglycosidic bond of the donor substrate is broken, the glycosidic bond to the acceptor Thr is formed, and a proton from Thr is transformed to the leaving oxygen of the phosphate group. Four atoms participate in the reaction: anomeric carbon at position 1 of the GalNAc ring (C1), oxygen of the phosphate group (O1), oxygen of the acceptor threonine (OA), and hydrogen of the acceptor threonine (H). The atoms are shown in Figure 1. To describe the reaction, we chose two sets of collective variables (CVs): one for the metadynamics (MTD) and one for the string method (STM). For the MTD, the computational cost grows exponentially with the number of collective variables; it is, therefore, crucial to use as few of them as possible. In our case, we utilized two distance differences as CVs. The difference between the distances C1−OA and C1− O1 was used to describe the breakage/formation of the glycosidic bond. The difference between O1−H distance and OA−H distance was used to describe the H transfer. The advantage of the string method is that the exploration of the reaction can be performed in a multidimensional space of collective variables. Therefore, for the STM simulation, we used the four distances (C1−OA, C1−O1, O1−H, OA−H) as four independent CVs. STM Details. An in-house implementation of the String Method in PMFlib43 was employed for the STM simulation. The initial path estimation was taken from our previous study using the nudged elastic band method. The same four distances (given above) were used for the path description, and the same number of beads (30) were utilized. The system was driven along the path from the equilibrated structure to generate the initial bead positions. The optimization step of STM consists of an equilibration stage and an accumulation stage. The equilibration stage consists of 200 MD steps, allowing the system to adapt to the updated bead positions. The accumulation stage is 2000 MD steps long, and, during this period, the forces acting on the bead in the direction perpendicular to the path are accumulated. Based on these forces, new bead positions are calculated, and the next round of optimization is performed. Between the optimization steps, the positions of the beads are allowed a maximum of 0.05 Å movement. Once the path optimization is completed, a production run of 25 000 MD steps (3 ps) is performed, and the final forces/free energy profile are calculated. Error estimation for the free energy profile was calculated by integrating the estimated standard uncertainties of the PMF forces on individual beads (including the effect of time correlations). The String Method internally describes the reaction path by a parametric spline curve with a parameter “alpha” going from zero (reactant) to one (product). To allow a direct comparison with the NEB-based path, we assigned an “alpha” value to each NEB image by projecting it onto the spline curve. The “alpha” value then corresponds to a point on the STM spline curve that has the lowest Euclidean distance in the CV space to a given

kcal/mol restraints were 50 and 70 ps long, respectively. The final stages of equilibration were a simulation with a 1 kcal/mol restraint placed on the protein backbone and a simulation with no restraints, both 300 ps long. A CUDA-enabled version of pmemd34 from the AMBER35 package was used to perform classical force field molecular dynamics on the GPU. The equilibrated structure was simulated for 30 ns at 300 K and 1 bar. Berendsen weak coupling thermostat36 and Berendsen barostat were used. QM Region. The system was divided in such a way that only atoms relevant to the reaction studied were treated at the QM level. The QM zone composition is described in Table I and Table I. Residues Included in QM Regions in NEB and CPMD Simulations included part residue

NEB

CPMD

function

Arg208 Asp224 His226 Ala307 Gly308 Gly309 Trp331 Glu334 His359 Arg362 His365 Tyr367

side chain side chain side chain without NH whole without CO side chain side chain side chain side chain side chain side chain

N/A side chain side chain N/A N/A N/A side chain N/A side chain side chain without CB N/A side chain

GalNAc H-bond25 Mn ligand21 Mn ligand21 GalNAc binding pocket21 GalNAc binding pocket21 GalNAc binding pocket21 GalNAc CH-π26 GalNAc H-bond21 Mn ligand21 phosphate H-bond23 W331 NH-π phosphate H-bond

shown schematically in Figure 1. The QM zone used in NEB, taken from our previous study,30 was made to include a large portion of the enzyme active site to account for other less critical interaction. The large zone was possible due to the lower computational cost of the calculations. Mg2+ ion was used in place of the native Mn2+. The expected catalytic activity with Mg2+ is 89% of the native activity.27 We have also previously verified that using Mg2+ does not significantly alter the results of QM/MM studies.30 Table I shows ppGalNAcT2 protein residues included in the QM zone and their roles; also, four water molecules were included in the smaller QM zone and an additional two in the larger one. The coordination sphere of Mg2+ cation consists of Asp224, His226, and His359 residues, one water molecule, and two oxygens of the diphosphate Due to the high computational cost, it is not possible to run QM/MM CPMD simulation with the full QM zone used in the NEB study. Therefore, only residues and water molecules either coordinating the metal cation or directly interacting with the phosphate part of the substrate were included. The hydrogen capping method was used to saturate the QM region in the CPMD. The total number of QM atoms used in CPMD (including the capping atoms) was 145 and 275 in the NEB calculations. CPMD Simulation. The simulations were conducted using CPMD version 3.13.237 with the Gromos QM/MM interface.38 The electronic system was described using the PBE39,40 functional with Grimme’s DFT-D2 dispersion correction41 and Troullier-Martins norm-conserving pseudopotential.42 The time step used was 5 au (0.12 fs) and fictitious electron mass 600 au. The plane-wave cutoff used was 70 Ry. The MM equilibrated system was subject to quick temperature quenching to alleviate any clashes between atoms when 6064


Article


Figure 1. Scheme of QM zones used in CPMD (black) simulations and NEB (gray) simulations. Arrows indicate the reaction.

previously used distance-difference scan coordinates, new combinations of scan coordinates were applied: • R−R′: Coordinate #1 = C1−OA distance, Coordinate #2 = difference between the O1−H and OPA−H distances • P−P′: Coordinate #1 = C1−OA distance, Coordinate #2 = difference between the O1−H and OD1−H distances Target values of these coordinates were extracted from the MTD structures representing R′/P′ shown in Figure 6. Harmonic restraints acting on these coordinates were gradually shifted by 0.3 Å after each geometry optimization run until the target values were reached. The restraints were then removed, and a fully unrestrained geometry optimization of the candidate R′/P′ was carried out with tight convergence criteria (maximum gradient element less than 0.001 hartree Å−1) to accurately locate the respective stationary point. Finally, a NEB optimization run between R and the newly located R′ (or P and P′) was carried out until the maximum element of the overall NEB gradient decreased below the default threshold of 0.05 eV Å−1 (0.0018 hartree Å−1). To test the ability of various density functional approximations to describe the glycosylation reaction correctly, single point energy evaluations were carried out on the main optimized path. The following density functionals were tested: BLYP,50−53 B3LYP,54 BP86,50,55,56 M06-2X,57,58 PBE,39,40 PBE0,59,60 S12g/S12h,61 and TPSS.62,63 Normal Mode Analysis. The vibrational analysis was conducted on the optimized reactant, product and transition

NEB image. This distance is plotted in Figure S1 to show that the projection is physically meaningful (All distances are low enough; thus, the original and projected points represent the same “chemical state”.). MTD Details. The strategy for the MTD simulation was to obtain a rough free energy landscape relatively quickly. Accordingly, the MTD strategy we employed is relatively crude. The width of the Gaussians deposited was 1.0 au (0.529 Å), the height was 0.004 au (2.51 kcal/mol), and the deposition time was 400 steps (48.4 fs). Wall boundaries were placed at values of 6 au (3.18 Å) for CV1 and 6.5 au (3.43 Å) for CV2 to keep the system from dissociating. The electronic degrees of freedom were thermostated based on an unbiased simulation. NEB Details. The original optimized NEB path published in our previous study30 was extended in both directions to investigate two new minima (R′, P′) discovered by MTD. The protocol used for the extensions is the same as for the main path, and we will thus reiterate it only briefly. All NEB calculations were done using ADF44−46 version 2013.01c coupled with a locally modified version of the Python ASE toolkit. The OPBE39,40,47 density functional with DFT-D3 dispersion correction48 and TZP basis set49 was used throughout the NEB study. First, initial guesses of the new NEB path segments were generated by successive restrained geometry optimizations starting from the optimized R and P structures. Because neither of the transitions (R′−R and P−P′) is described well by the 6065


Article


which was performed in full Cartesian space, the String Method was performed in the space of the four CVs described in Methods. Based on the forces acting on the constraints, the position of the beads are obtained in each step until a final converged path is obtained. The optimized potential energy path served as an initial guess for the String Method. The optimization was run for 30 steps before the path was considered converged. The change in the CV values through the path optimization process is shown in Figure S2. It can be seen that CV1 and CV2, describing the breakage of the original glycosidic bond and formation of the new one respectively, converge quickly to their final values. CV3 and CV4 describing the proton transfer have more trouble converging; they were essentially jumping between two possible values at the start of the path from one optimization step to another. Figure S3 shows the preliminary free energy profiles obtained at each optimization step. The initial path guesses taken from the NEB study show a profile with 27 kcal/mol free energy barrier and 12 kcal/mol reaction free energy. During the optimization, the free energy barrier decreases below 20 kcal/ mol, and the reaction free energy reaches 0 kcal/mol. The accumulation phase of the STM optimization is only 240 fs long resulting in a limited sampling of the free energy and presumably rather big errors. We, unfortunately, do not have error estimations of the optimization steps. Figure S3, however, shows the standard uncertainties of the final production profile. Most of the preliminary profiles from optimization steps fit within them. The errors on the preliminary profiles are expected to be bigger. As such the free energy profiles of the optimization steps are not a reliable way to judge convergence. More important to us are the changes in the geometry of the path. This can be judged by looking at the progress of the CVs during optimization shown in Figure S2. Based on the CVs changes, we decided to terminate the path optimization after 30 steps. It seems unlikely that further optimization would improve the path dramatically and, in our opinion, would not be worth the computational time. The final reaction profile obtained from the production run on the final path is shown in black and gives a 22.5 ± 7.6 kcal/ mol free energy barrier and 5.2 ± 10.0 kcal/mol reaction free energy (95% confidence intervals). It is worth mentioning that the integration error accumulations along the profile mean that the reaction free energy is more burdened, decreasing its reliability. The final free energy profile including integration error estimation is shown in Figure 2, together with the

state stationary points from our previous study. Hessian matrices were generated by numeric differentiation in all Cartesian coordinates, using two-point central finite differences of the analytic gradients with a differentiation step of 0.02 Å. These gradient evaluations were done using ADF 2013.01. The calculated matrix was symmetrized by averaging it with its transpose. Normal mode vibrations and thermodynamic contributions were then calculated using the “aoforce” and “freeh” programs from TURBOMOLE 7.00. Coupled Clusters Calculations. All calculations were performed using ORCA64 version 3.0.3. Single-point energy calculations were carried out on the 30 NEB images of the neutral 275-atom QM region from our previous QM/MM study.30 No periodic boundary conditions or continuum solvation models were applied. Spin-restricted Hartree−Fock energies were calculated using the def2-SVP and def2-TZVPP basis sets65,66 in combination with the Split-RI-J charge fitting67 and Chain-Of-Spheres (COS) exchange fitting methods.68 The “very tight” preset for SCF convergence criteria was applied to prevent numerical noise from spoiling the results. Optimized Kohn−Sham orbitals from a preliminary calculation using the PBE density functional39,40 and def2-SVP basis set were used as the initial orbital guess for HF/def2-SVP. Initial guesses for calculations using larger basis sets were obtained by projection of the HF orbitals optimized using the closest smaller basis. Second-order Møller−Plesset correlation contributions were subsequently calculated using the optimized RHF/def2-SVP and RHF/def2TZVPP wave functions. The highest level of correlated calculations was represented by CCSD(T) coupled cluster theory (including all single and double excitations and a perturbative correction for triple excitations). The reference wave functions used were the same as in the MP2 calculations. To make these computationally demanding calculations tractable on such a large system, we chose the DLPNO18 (Domain-based Local Pair Natural Orbitals) linear-scaling approximation. All the screening and cutoff parameters of DLPNO were kept at their conservative default values. Complete basis set extrapolations were carried out based on the def2-SVP and def2-TZVPP energies using two schemes. In both cases, the extrapolated correlation contribution was added to an RHF/CBS energy calculated using the established exponential formula69 • A direct extrapolation of the DLPNO-CCSD(T) energies was carried out using the power scheme70,71 • Alternatively, the MP2 correlation energy was extrapolated to the CBS limit using the power scheme, and a single-point coupled clusters correction was added Optimized values of the alpha and beta parameters for the pair of basis sets used were taken from ref 72.

■

RESULTS AND DISCUSSION STM Results. In our previous work30 we utilized the Nudged Elastic Band methodology to find the minimum potential energy path. We also obtained the optimized structure of the transition state of the ppGalNAcT2 glycosyltransferase reaction. The aim of this work is to extend the study into the state of the free energy space. The String Method is a path optimization method similar to NEB. The reaction in String Method is represented by so-called beads that are connected to form the path. Each bead is a separate MD simulation constrained in the CV space. In contrast to the NEB method,

Figure 2. Final free energy profile from STM (green) and optimized profile from NEB (purple). Standard uncertainties of the STM free energy profile are shown as error bars. 6066


Article


Figure 3. Evolution of relevant distances along the reaction path. (A) Formation of a new C1−OA glycosidic bond (green − STM, black − NEB) and breakage of an old C1−O1 glycosidic bond (blue − STM, purple − NEB) and (B) H transfer from OA−H (brown − STM, yellow − NEB) to O1−H (orange − STM, red − NEB).

potential energy profile obtained previously. The force acting on the STM beads is shown in Figure S4. The force crosses zero around the bead corresponding to the highest point on the free energy profile. We consider this bead to represent the transition state of the reaction. The shape of the free energy profile is comparable with the profile from NEB. Perhaps most importantly, the positions of the maxima corresponding to the reaction barrier are the same. However, the free energy profile is energetically higher and less flat. The first shallow minimum present in the NEB profile after the initial barrier is not observed in the free energy profile. Instead, the energy steadily rises until it reaches the real reaction barrier. The reaction barrier is also higher, 22.5 ± 7.6 kcal/mol from the STM compared to 14.1 kcal/mol from NEB. Comparing the free energy profile from STM with a rigid rotor−harmonic oscillator−ideal gas estimation of free energy on the NEB stationary points in Figure S5 shows that the magnitude of the entropic contribution to the barrier height is about 4.5 kcal/mol. However, the estimation indicates that the barrier would be lowered, contrary to the higher barrier observed from STM. The evolution of the STM CVs along the path (shown in Figure 3) gives us a more detailed understanding of the reaction. We can see that the reaction proceeds by breaking the original glycosidic bond and gradually extending the C1−O1 distance, while the C1−OA distance decreases at a slower pace. The reaction, therefore, goes through an oxocarbenium ion pair state. The structure corresponding to the highest point on the free energy profile is shown in Figure 4. The C1−OA distance decreases to the point of forming the new glycosidic bond. The distances corresponding to the H-transfer see relatively little change until the transfer itself occurs simultaneously with the formation of the new Thr-GalNAc glycosidic bond. It is apparent from Figure 3 that the distances along the free energy path differ only slightly from the ones obtained along the potential energy path, and the nature of the reaction is the same. The biggest difference is in the C1−O1 distance of the dissociating glycosidic bond around alpha = 0.7. In the NEB path, this separation reaches 3.6 Å and then settles at 3.4 Å in the product state, while the C1−O1 distance goes smoothly to 3.4 Å and maintains this distance to the product state without any bump in the STM path. This part of the profile corresponds to the transition state. The comparison of the transition state structures between the NEB and the STM is of particular interest. Figure 4 shows the comparison of the optimized TS from our previous work with an average structure corresponding to the STM bead with the highest free energy.

Figure 4. Comparison of optimized NEB Transition State structure (green) and the ensemble average structure corresponding to the highest energy bead along the STM path (color). Hydrogen bonds are highlighted in red and blue.

The distances are listed in Table II. The TS structures are quite similar. The differences in the relevant distances are 0.05 Å or below, with the sole exception being the C1−O1 distance with the difference of 0.19 Å discussed above. The source of this difference is not entirely clear but may be a result of entropic effects not accounted for in the static TS optimization. Other geometrical differences are in the orientation of the N-Acetyl group and the OH6 group of the GalNAc. The orientation of the OH6 group is of little interest to us. The N-Acetyl group orientation results in a difference in hydrogen bonding to the leaving phosphate (shown in Figure 4). This may be linked to the difference in the distance of the GalNAc ring to the leaving phosphate between the TS structures, which appears in the C1−O1 distance comparison. The consequences of this difference for drug design purposes are unclear. MTD Results. A metadynamics simulation was performed to generate a rough free energy landscape of the reaction. It can serve to verify that the optimized reaction path lies in an energetically favorable region and represents a realistic path from reactants to products. Due to the dependence of the computational demand on the number of CVs, the simulation was performed in a 2D space using the distance difference CVs described in Methods. The distance CVs used in the STM can be easily transformed into distance difference form, and the 6067


Article


Table II. Relevant Distances and Cremer-Pople Ring Puckering Parameters in NEB and STM Reactant (R), Transition State (TS), and Product (P) Structuresa R C1−OA (Å) C1−O1 (Å) O1−H (Å) OA−H (Å) Cremer-Pople ϕ (deg) Cremer-Pople θ (deg) Cremer-Pople R a

TS

P

NEB

STM

NEB

STM

NEB

STM

3.05 1.51 1.78 0.97 246.8 9.1 0.596

3.14 1.50 1.83 0.96 205.9 7.2 0.567

2.33 3.60 1.40 1.04 249.8 44.4 0.558

2.35 3.41 1.45 1.04 253.9 44.9 0.537

1.47 3.48 1.04 1.43 259.0 24.3 0.553

1.44 3.52 1.03 1.50 235.9 19.2 0.549

The STM distances are the final optimized values; the STM Cremer-Pople parameters are ensemble averages from the corresponding bead.

this minimum is thus probably overestimated. The other noteworthy point in the reactant space is marked R and represents the starting point of our previous NEB study and this STM simulation. This point corresponds to a state in which the threonine acceptor forms a hydrogen bond with the leaving oxygen O1 of the phosphate group, and the C1−OA distance is shortened to ∼3.1 Å. The free energy at the R point is ∼−48.8 kcal/mol. From the R point, the reaction proceeds through a mostly horizontal valley connecting the reactant and product minima. The valley stretches from ∼1.8 Å to ∼−2.5 Å in the space of CV1, representing the breaking/formation of the glycosidic bond, and then from ∼1 Å to ∼−1 Å in the space of CV2, describing the H transfer. The H transfer itself occurs simultaneously with the formation of the new glycosidic bond near the product state P. The point of highest free energy across the reaction path (labeled TS′) has a value of −30.4 kcal/mol. The product state is in the bottom half of the surface and is characterized by the presence of two minima. The minimum labeled P corresponds to the product in both the NEB and STM simulations. The value of the free energy at the P minimum is ∼−54.0 kcal/mol. The additional minimum P′ differs mainly in the space of CV2 and corresponds to the different orientation of the now protonated leaving phosphate group. The minimum labeled P represents a conformation in which the protonated O1 atom forms a hydrogen bond with the OA atom. The minimum labeled P′ represents a conformation in which the protonated phosphate group is pointed toward the Mg2+-coordinating oxygen of Asp224. This minimum is not relevant and may be the result of artifacts caused by the usage of distance difference as a CV. Details are provided in the Supporting Information. The free energy barrier of the reaction is ∼32.0 kcal/mol when we consider the point R′ as a reactant. R′ may, however, be too oversampled and artificially too deep due to the presence of wall boundaries, as we discussed above. When we only consider the part of the free energy surface corresponding to the paths in the STM and NEB simulations and look at the reaction from R to P, the barrier is ∼18.4 kcal/mol. It should be noted that the reconstructed free energy surface does not represent a fully converged free energy landscape, as that was not our main goal. The reaction barrier was crossed only twice, and the MTD setting used was rather crude. The error in the free energy should be at least equal to the height of the Gaussians used, i.e. 2.5 kcal/mol (probably more). Even then the obtained free energy barriers correspond to the values expected from the enzymatic reaction and are in line with the experimental rate constant20 and previous computational

path can, therefore, be easily mapped onto the surface. The production metadynamics simulation was run for a total of 18 ps, which was the time it took to cross the reaction barrier, explore the conformational space of the products, and cross back to the reactants space. At that point, the simulation was stopped. During the simulation, 366 Gaussians were deposited. The reconstructed free energy surface is shown in Figure 5, and the evolution of the two CVs is shown in Figure S6. Points

Figure 5. Reconstructed free energy surface from metadynamics. The optimized NEB and STM paths are shown in blue and black, respectively. The extension of the original NEB path is shown in lighter blue.

of interest are marked in Figure 5, and their representative structures are shown in Figure 6. The reconstructed surface features a very broad minimum corresponding to the reactant state of the reaction. The wide nature of the minima can be in part attributed to the intrinsic flexibility of the distance differences chosen as CVs. The lowest point (marked R′) represents the state in which the threonine acceptor develops a hydrogen bond with OPA of the phosphate group and with the C1 and OA atoms, separated by ∼4.1 Å. This point corresponds to the lowest free energy point at −62.4 kcal/mol. It should be noted that this point lies in a corner formed by the wall boundaries keeping the system from dissociating. The depth of 6068


Article


Figure 6. Representative structures showing the progress of the reaction. The structures correspond to the labeled point in Figure 5. The R, TS, and P structures are taken from STM, and R′, TS′, and P′ are taken from MTD. Shown is just a subset of the QM zone including acceptor Threonine (ThrAcc), GalNAc, uridine diphosphate (UDP), and Mg2+ cation; where relevant also Asp224 coordinating Mg2+ and Trp331 stabilizing the leaving phosphate group is shown.

results. Barrier heights and reaction energies calculated using various methods are compared in Table III. The main purpose of the MTD was to serve as a qualitative description of the glycosyltransferase reaction and provide context and verification for our STM study. STM gives a better and more detailed description of the optimized reaction path

but cannot explore the whole free energy landscape of the reaction. The free energy surface from MTD, although rough, can determine whether the path we found truly lies in the global minimum or became trapped in a local minimum path. The optimized STM path we obtained follows the free energy valley reconstructed from MTD simulation, as can be seen in 6069


Article


and OA atoms gradually shortens (structure TS in Figure 6). Finally, the new glycosidic bond with threonine is formed, and the proton is simultaneously transferred from the threonine to the oxygen of the phosphate group (structure P in Figure 6). The now protonated phosphate group can reorient itself toward the Asp224 involved in coordinating the Mg2+ ion (structure P′). Further stabilization is provided by the N-acetyl group of the donor substrate and Trp331. Ring Puckering. An important feature when discussing glycosyltransferase reactions is the itinerary of ring puckering conformations of the reaction. The previous study gave us insight into the conformational behavior of the GalNAc ring in the context of NEB path optimization. For the STM simulation, we obtain an ensemble of structures for each point on the free energy profile. These results can complement the information from static NEB structures by assessing the dynamic behavior and flexibility of the ring puckering along the reaction path. Figure 7 shows the itinerary of the Cremer-Pople ring puckering parameters73 from the optimized NEB structures (in

Table III. Comparison of Energies of the (Estimated) Transition and Product States Relative to the Reactant Statea method ΔEpot

ΔG

a

NEB NEB NEB Gómez et al.28 (1D scan) NEB STM MTD Lira-Navarrete et al.29 (MTD)

level of theory

transition state

product

OPBE M06-2X DLPNOCCSD(T) M05-2X

13.8 14.1 13.2

−6.7 −11.3 −10.9

19.8

0.3

OPBE PBE PBE PBE

8.1 22.5 ± 7.6 ∼18.4 ∼13

−7.4 5.2 ± 10.0 ∼− 5.2

All values are in kcal/mol.

Figure 5. The points and structures labeled R, TS, and P correspond to the starting point, highest energy point, and final point of the optimized path, respectively. This agreement between the STM and MTD indicates that the path optimization indeed found the correct reaction path. The slight difference between the STM and NEB paths can also be seen from the projection of the path onto the two-dimensional distance difference space of the MTD simulation. The NEB path travels further along the horizontal axis describing the glycosidic bond before dropping to the product state. This corresponds to the different behavior of the C1−O3B distance discussed above. To check the relevance of R′ and P′ minima identified in MTD simulation, the original NEB was extended in both directions (shown in Figure 5 as lighter blue). The resulting potential energy profile is shown in Figure S7. The extension toward R′ shows a slight rise in the potential energy instead of expected descent toward a deep minimum. Additionally, even though the R′ end of this path was initially placed at the position of the corresponding minimum on the MTD FES surface, it moved considerably closer to R during the subsequent geometry optimization. Still, the hydrogen bonding arrangement of R′ was fully retained. We assume this as an indication that the R′ minimum seen in the MTD simulation is, in fact, an artifact of the simulation. Extending the NEB path toward P′ shows another minimum slightly higher in energy than P. There is, however, a large potential energy barrier separating the two minima. The presence of a barrier comparable to or even higher than the barrier of the enzymatic reaction itself questions the relevance of such a minimum. We attribute the P′ minimum seen in MTD to an artifact caused by the use of the distance difference collective variable (further discussed in Supporting Information Figures S8 and S9). Reaction Description. Figure 6 illustrates the reaction process and the interactions involved. The reaction proceeds by breaking the original glycosidic and separating the C1 and O1 atoms to about ∼3.5 Å, while keeping the C1 and OA atoms at a moderate distance of ∼3.2 Å (structure TS′ in Figure 6). The leaving O1 forms a hydrogen bond with the acceptor threonine, but the proton transfer does not occur at this point. Further stabilization of the leaving phosphate group is provided by the N-acetyl group of the donor substrate, Trp331, and water molecules, all via hydrogen bonding. The separation of the C1 and O1 atoms is maintained, while the distance between the C1

Figure 7. Evolution of GalNAc ring puckering along the NEB (red) path and the STM (blue) path.

red), the production run ensemble averages corresponding to the beads along the STM path (in blue). The standard deviations of the puckering angles are shown in the Supporting Information (Figure S10). Both NEB and STM start from the same region around the 4C1 chair conformation. The GalNAc ring exhibits a large degree of flexibility in the reactant state. This flexibility decreases significantly as the glycosidic bond starts elongating. Interestingly, the ring conformation changes differ between the STM and NEB pathways. The optimized NEB structures follow straight from the 4C1 to the 4E conformation, while in the STM simulation the ring conformation changes along the reaction are more disordered and reach a shape close to the 4H3 conformation. Importantly, the transition state structures (labeled TS in Figure 7) exhibit the same ring conformation in both the NEB and STM studies. The low standard deviation in the transition state ensemble 6070


Article

Journal of Chemical Theory and Computation (Figure S10) revealed a quite rigid ring conformation and thus indicates an important feature of the glycosyltransferase reaction. This suggests that mimicking the proper ring shape is an important structural feature in designing potent transition state analog inhibitors for glycosyltransferases. From the transition state, both the NEB and STM follow roughly along the same ring conformation to a more flexible product state with a ring conformation slightly different from the reactant state. CBS Extrapolation of DLPNO-CCSD(T) Energies. Correlated wave function theory methods in general exhibit a relatively slow convergence of the calculated energies toward the complete basis set limit. There are several widely accepted extrapolation schemes solving this issue for pure coupled cluster calculations. However, their applicability to the results calculated using the DLPNO approximation has not been verified yet. Figure S11 shows the part of the reaction energy profile corresponding to the dissociation of the donor glycosidic bond. This covers the transition from an equilibrium geometry to a dissociated oxocarbenium-like state. Because of the corresponding large change in the correlated pairs, this process should offer a real test of the reliability of the DLPNO approximation for the studied system. Comparison of the calculated energy profiles reveals that using the DLPNO-CCSD(T) method does not introduce any significant artifacts compared to the established MP2/CBS extrapolation. This holds even if the DLPNO-CCSD(T) correlation energy is extrapolated directly instead of merely using it as a correction to MP2/CBS correlation energy. As the MP2/CBS calculations are much more computationally demanding than direct DLPNO-CCSD(T)/def2-TZVPP calculations, it seems reasonable to avoid doing MP2 completely. We have thus opted for the direct extrapolation approach in the rest of this study. Of course, one may wonder if using solely double- and tripleζ basis sets is enough to obtain reasonable CBS extrapolation results. Unfortunately, running triple-ζ MP2 or DLPNO-CC calculations is the best we can do with available computational resources for such a large system. Quadruple-ζ HF would be tractable but would not produce any improvement, since the TZVPP energy profile is already indistinguishable from the CBS-extrapolated one (see Figure S12). We can thus only offer a crude estimate of the correlation energy lost due to basis set incompleteness. By comparing the correlation energy profiles in Figure S13, we can estimate the worst-case correlation energy extrapolation error to around 1 kcal/mol. Even such an error would be acceptable for our study, since we are comparing the energies to DFT-based results that are certainly no more accurate. Additionally, there are other important sources of correlation energy errors (for example higher-order coupled cluster excitations) that we also cannot eliminate due to resource constraints. Comparison of CC and DFT Energy Profiles. Figure 8 shows a comparison of the final CBS-extrapolated potential energy profile with the results of our previous DFT study. It is reassuring that the DFT results compare very favorably with the new coupled cluster reference data. The M06-2X meta-hybrid functional,57,58 in particular, exhibits an almost unbelievable accuracy in describing the heights of both reaction barriers and the shallow local minimum corresponding to the oxocarbenium intermediate. This success is probably due to the heavily parametrized nature of this Minnesota functional.

Figure 8. Comparison of QM energy profiles along the NEB reaction path.

On the other hand, the popular B3LYP hybrid density functional54 consistently underestimates the barrier heights by about 2 kcal/mol, although it still exhibits reasonable qualitative agreement with the reference. However, this only applies when the DFT-D3 dispersion correction is not used, otherwise the agreement would be considerably worse. In contrast, the OPBE density functional39,40,47 combined with DFT-D3 gives a slightly worse qualitative description of the reaction, where the positions of barriers are somewhat shifted, and the intermediate is only distinguishable on a separate plot of QM energy (Figure S14). However, the overall shape of the potential energy profile is still preserved, and the barrier height is reproduced with an error of less than 2 kcal/mol. A variety of other density functionals were also tested, and the results are shown in Figure S14. It is somewhat surprising that employing the DFT-D3 correction does often lead to worse results than with an uncorrected functional. This is particularly true of most pure GGA functionals, with OPBE being the sole exception. Even the S12g functional61 that is closely related to OPBE and designed to be used with dispersion correction performs better when DFT-D3 is switched off. We have also tested the performance of two GGA functionals with and without DFT-D3 in ORCA to make sure that the disappointing performance of GGA is not due to an implementation problem in ADF (Figure S15). However, the results from ORCA are qualitatively the same, with differences around 2 kcal/mol that can be attributed to the different approximations, basis sets (STO vs GTO), and numerical integration methods used in both software packages. Timings. Computational costs associated with methods used in our study are shown in Table IV regarding CPU hours (assuming common midrange Sandy Bridge Xeon CPUs like for example E5-2670). The most costly method is STM with a total of 690 000 CPU hours needed to optimize the pathway and obtain final free energy profile. The parallelization of the STM between in our case 30 independent beads makes this method usable provided there be access to supercomputing facility. The metadynamics has much lower total computational cost, but it was not parallelized in the way STM was. There is an option of using multiple walker approaches to speed up the MTD simulation. However, it has not been tested in the context of QM/MM AIMD simulations. The NEB calculation provides a very cost-effective alternative. Just as STM, NEB is also trivially parallel because the gradient evaluations on individual NEB images (beads) are independent. Even with a 6071


Article


various approximations used in CPMD and ADF. Even the DFT-D3 dispersion correction may be playing an important role. Its contribution to the NEB profile using the PBE functional lowers both the barrier height and reaction energy by 4−5 kcal/mol (Figure S14). It is, however, unclear how well the older DFT-D2 correction used in the CPMD STM simulation matches these results. Perhaps more important is the geometrical comparison between the methods. The relevant distances follow the same evolution along the path, apart from the distance between the C1 carbon and dissociating oxygen of the phosphate group around the maxima of the energy profile. This difference manifests itself in the distance between the paths (Figure S1) and in the structure of the transition state. The optimized transition state structure and the ensemble average structure of the transition state from STM are otherwise the same. There are multiple possible sources for this difference. One explanation would be the effect of entropy, which is explicitly included in the STM. Another possibility is the effect of the functional. Another interesting feature to compare is the conformational behavior of the GalNAc ring along the reaction path. Both static structures from NEB and ensemble averages from STM agree in the structures of the transition state and product. Their conformational behavior preceding the TS structures is, however, different. In NEB, the optimized structures along the path transition proceed from the initial 4 C1 toward the 4E conformation and then to the transition state. In the STM, the GalNAc ring prefers structures closer to 4 H3 conformations and skips over the 4E part of the conformational space. The structures obtained from NEB represent a minimum in the potential energy space. In STM, each bead is restrained in the CV space to the minimum free energy path, but the rest of the system not included in the CVs is free to move. This results in a difference in surroundings of the GalNAc moiety, that may favor ring conformations away from the potential energy minima seen in NEB. An important result obtained from STM is the ensemble average, as well as the standard deviation of the ring puckering parameters at the transition state, shown in Figure S10. From the standard deviation, we can see the increased rigidity of the ring conformation at the transition state. Maintaining the ring conformation may, therefore, be a major consideration in the design of transition-state mimetics for rational drug design. We performed metadynamics simulations to obtain a free energy surface of the reaction. The MTD complements the path optimization methods by providing an overview of the whole reaction space. This enables us to check the validity of the paths obtained. Path optimization methods carry a risk that the reactant, product, and/or transition states get stuck in a local minimum. The results then describe an incorrect path and may, therefore, misrepresent the energetic barrier and the structure of the transition state. The reaction paths we obtained from both NEB and STM nicely follow the minimal free energy region connecting the reactant and product minima. The product state P, the end point of both paths, sits in the corresponding free energy minimum. The highest point on the free energy surface along the reaction valley TS′ is near the transition state TS from STM. However, this region of the surface is less sampled, as the system only crossed the barrier twice. The size of Gaussians used also limits the resolution of the resulting surface, making precise localization of the transition state difficult. One advantage of STM is that every point along the reaction path, including the transition state, is

Table IV. CPU Timings of the Methods Used method

stage

MTD STM

N/A path optimization production run path optimization single-point DLPNOCCSD(T)

STM NEB NEB NEB

total CPU hours 40 000 500 000

details

190 000

30 beads × 30 optimization steps × 2200 CPMD steps 30 beads × 25000 CPMD steps

9 400

28 active beads × 100 iterations

3 100 19 200

30 beads, M06-2X 30 beads (CPUh breakdown: 50 HF/ CBS + 80 CC/SVP + 500 CC/ TZVPP; optionally 115 MP2/SVP + 1405 MP2/TZVPP)

large QM region, tens of NEB optimization iterations can be done in a single day on relatively small computational clusters. The high-quality DLPNO-CCSD(T)/CBS calculations are also within reach of many computational groups even without access to large supercomputers. However, the timings shown in Table III were done on an SGI UV2000 machine with 20 GiB RAM per one CPU core. Although the calculation can probably be successfully done even on computers with much more limited RAM resources, the CPU time needed may then increase considerably.

■

DISCUSSION We compare the results from a static path optimization method (NEB), with results from a free energy path optimization (STM). Both methods agree to the general shape of the energy profile and the position of the maxima along the path. Relative energies of the stationary points are compared in Table III. The barrier height from STM is 22.5 ± 7.6 kcal/mol (95% confidence interval), which is higher than the 13.85 kcal/mol barrier from NEB. The value of 22.5 ± 7.6 kcal/mol lies in the expected range, given the experimental estimation.20 Where the NEB and STM profiles disagree is the product state. In the potential energy profile, the product lies lower than the reactant, with a reaction energy of −6.70 kcal/mol. The STM, on the other hand, has a product higher in energy and a reaction free energy of 5.2 ± 10.0 kcal/mol. Positive values of reaction free energy from STM simulation are something we also observed in another system.16 A partial source of such a difference could be the integrated sampling error in the STM free energy profile, which increases along the path. The product is therefore more burdened by the error than the reactant. The cause of this error lies in an insufficient sampling of the PMF force. Increasing the simulation time of the production run beyond 3 ps/bead would, however, lead to unreasonably high computational costs. We have also attempted to explain the role of entropic contributions in the observed differences between STM and NEB by estimating these contributions from normal mode theory. Unfortunately, the results show that the free energy barrier at the optimized NEB-based transition state would be about 4.5 kcal/mol lower than the potential energy barrier. The most probable explanation for this discrepancy is that the harmonic oscillator approximation does not hold for the many weak nonbonded interactions and hindered rotations present in this biomolecular system. A significant part of the quantitative differences between the NEB and STM energies may also be due to the different QM implementations, basis sets, QM/MM coupling schemes, and 6072


Article


much more crucial in methods such as metadynamics, where one is usually limited to just 2−3 CVs. Coincidentally in our case, the four distances describing the reaction can be transformed into two distance differences making it suitable for MTD simulation. One can go even further and combine three distances into one CV as demonstrated by Gómez et al. in the study of the same glycosyltransferase.28 However, even simple distance difference CVs such as the ones we used can lead to artificial behavior, demonstrated in the Supporting Information (Figures S8 and S9). Another option is to use CVs based on coordination numbers, the approach used by LiraNavarrete et al. to once again study ppGalNAcT2.29 The coordination numbers make the comparison with the distance based STM difficult and were not suitable for our study. The relevance of our results regarding the MTD simulation to other systems is, therefore, limited. However, it should apply to other retaining glycosyltransferases, other reactions that proceed through an SNi-like mechanism, or, to a lesser extent, to any enzymatic reaction involving a simple transfer between two substrates coupled with proton transfer. When we compare the QM description offered by approximate density functionals with reference DLPNOCCSD(T)/CBS energies, we can conclude that DFT offers a reasonable level of accuracy. In particular, the M06-2X metahybrid functional57,58 matches the coupled clusters reference profile almost perfectly, with errors under 1 kcal/mol for all the stationary points. We can thus confidently recommend M06-2X as a standard for testing other density functionals on this class of reactions. However, such a computationally demanding density functional cannot be routinely used for molecular dynamics or NEB path optimizations. Surprisingly, even a pureGGA density functional like OPBE can describe the catalytic reaction of ppGalNAcT2 with sufficient accuracy. OPBE thus offers a very interesting choice for explorative studies of reaction mechanisms, with a computational cost orders of magnitude lower than the other methods discussed here. However, employing the OPBE functional for other glycosyltransferases should be done with caution.

sampled equally. The transition state and the free energy barrier can, therefore, be obtained more reliably. We can also gain more insight into the actual structure of the transition state and its dynamic behavior. An example of such information would be the ring puckering discussed above. There is, however, the slight discrepancy between the free energy surface from MTD and the paths from NEB and STM in the reactant state. The starting point of both paths lies on the edge of a broad minimum in the surface and not at the lowest point R′. Due to the presence of wall boundaries, the point R′ is oversampled and therefore artificially deeper. The wall boundaries are a necessary evil to stop the dissociation of the substrates so that we do not waste resources simulating a process we are not investigating. Additionally, the difference in hydrogen bonding between R and R′ is not adequately described by the distance difference coordinates used in MTD. The apparent minimum in R′ is thus probably composed of a mixture of similar partially dissociated states. The artificial nature of R′ minimum is supported by the potential energy profile of the extended NEB path (Figure S7) and by the fact that the optimized R′ stationary point lies in a considerably different position than expected from MTD (Figure 5). The issue remains of choosing the reactant point in the path optimization. A notable shortcoming of our MTD simulation is the short simulation time of just 18 ps and a rather crude Gaussian deposition strategy. It should be noted that during that period the system crossed the reaction barrier twice, from reactant to product and back. This fulfilled our aim for the simulation to explore the free energy surface of the reaction and confirm the validity of NEB and STM pathways. Of course, to obtain reliable free energy information on the reaction including the barrier, it would be necessary to perform a “proper” fully converged MTD simulation with carefully selected Gaussian deposition strategy until the diffusive motion in the CVs space is achieved. That is however not currently feasible in the realm of QM/MM AIMD simulations. One can perform either MTD simulation with crude Gaussian deposition strategy to more quickly recross the barrier and obtain a global overview of the reaction space at the cost of accuracy of the free energy information or, alternatively, more gentle MTD simulation with more accurate free energy information at the cost of reduced sampling across the barrier. Since the latter approach has already been applied to the ppGalNAcT2 reaction using QM/ MM CPMD method by others,29 we decided not to devote our efforts in this direction. In any future studies based on free energy, the recommended procedure would be to perform MTD first and use the rough free energy surface to select the end points for STM path optimization and to construct an initial guess of the reaction path. Alternatively, if computational cost is a decisive factor, one may obtain most of the same key insights by using just a potential energy-based NEB path optimization. The selection of suitable collective variables is an important choice one has to make when studying enzymatic reactions. The increase in a number of collective variables leads to exponential increase in computational time. This is however not the case for the string method. The computational cost of STM is independent of the number of CVs used. The reaction mechanism of ppGalNAcT2 is simple enough that we were able to describe it by just four distances. However, the method has been used to study more complex mechanisms described using more CVs.16 Our conclusion regarding STM should be transferable to other enzymatic systems. The choice of CVs is

■

CONCLUSIONS Looking at the description of the enzymatic reaction of ppGalNAcT2 by various computational methods, we can see that all of the approaches used reach the same general conclusions. In particular, the minimum energy reaction paths based on either potential or free energy provide a very similar qualitative picture, predicting almost the same location of the transition state. However, any path optimization method should be combined with a different approach providing an overview of the shape of the potential/free energy landscape. This is the only way to validate the physical relevance of the optimized path. In our study, we have shown that metadynamics can be used for this purpose, providing a rough estimate of the regions of interest. Unfortunately, there are many pitfalls of metadynamics that make it not reliable enough to be used as the only tool for studies of reaction mechanisms. For example, reactant and product basins can be oversampled due to the boundary conditions or the inadequacy of used collective variables, and transition state regions can be poorly refined due to a limited number of barrier crossings. On the other hand, the results from metadynamics are good enough to be used as suitable starting points for path optimization. 6073



■

Comparing the results of a free energy-based (STM) reaction path with its potential energy counterpart (NEB), we can conclude that entropy does play a visible role in the studied reaction mechanism. However, its influence is mostly limited to regions less attractive from a drug design standpoint. Although the overall height of the energy barrier does differ significantly between the free and potential energy, the key geometrical information about the transition state is essentially the same. We can thus confidently recommend simple potential energybased NEB path optimization as a suitable alternative to an extensive STM free energy study, as the amount of information lost is greatly outweighed by the computational time savings. The final reassuring conclusion provided by our study concerns the widespread usage of DFT for investigating the reaction mechanisms of glycosyltransferases. By comparing the results of density functionals of various levels of complexity with high-quality reference DLPNO-CCSD(T)/CBS data, we have found that even some of the computationally cheapest GGA methods provide a sufficient level of accuracy. At the same time, highly complex meta-hybrid density functionals such as M06-2X can match the reference data almost perfectly. They can thus be seen as an inexpensive path toward high-quality energy profiles.

■

REFERENCES

(1) Hart, G. W.; Slawson, C.; Ramirez-Correa, G.; Lagerlof, O. Cross Talk Between O-GlcNAcylation and Phosphorylation: Roles in Signaling, Transcription, and Chronic Disease. Annu. Rev. Biochem. 2011, 80 (1), 825−858. (2) Gill, D. J.; Clausen, H.; Bard, F. Location, Location, Location: New Insights into O-GalNAc Protein Glycosylation. Trends Cell Biol. 2011, 21 (3), 149−158. (3) Dennis, J. Glycoprotein Glycosylation and Cancer Progression. Biochim. Biophys. Acta, Gen. Subj. 1999, 1473 (1), 21−34. (4) Maddry, J. A.; Suling, W. J.; Reynolds, R. C. Glycosyltransferases as Targets for Inhibition of Cell Wall Synthesis in M. Tuberculosis and M. Avium. Res. Microbiol. 1996, 147 (1−2), 106−121. (5) Berg, S.; Kaur, D.; Jackson, M.; Brennan, P. J. The Glycosyltransferases of Mycobacterium Tuberculosis–Roles in the Synthesis of Arabinogalactan, Lipoarabinomannan, and Other Glycoconjugates. Glycobiology 2007, 17 (6), 35R−56R. (6) Belyi, Y.; Aktories, K. Bacterial Toxin and Effector Glycosyltransferases. Biochim. Biophys. Acta, Gen. Subj. 2010, 1800 (2), 134− 143. (7) Tvaroška, I. Atomistic Insight into the Catalytic Mechanism of Glycosyltransferases by Combined Quantum Mechanics/Molecular Mechanics (QM/MM) Methods. Carbohydr. Res. 2015, 403, 38−47. (8) Henkelman, G.; Jónsson, H. Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 2000, 113 (22), 9978−9985. (9) E, W.; Ren, W.; Vanden-Eijnden, E. String Method for the Study of Rare Events. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 66 (5), 52301. (10) 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 (4−6), 429−436. (11) Sattelle, B. M.; Almond, A. Less Is More When Simulating Unsulfated Glycosaminoglycan 3D-Structure: Comparison of GLYCAM06/TIP3P, PM3-CARB1/TIP3P, and SCC-DFTB-D/TIP3P Predictions with Experiment. J. Comput. Chem. 2010, 31 (16), 2932−2947. (12) 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 (51), 17142−17154. (13) 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 (10), 4694−4707. (14) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99 (20), 12562−12566. (15) Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. String Method in Collective Variables: Minimum Free Energy Paths and Isocommittor Surfaces. J. Chem. Phys. 2006, 125 (2), 024106. (16) Kumari, M.; Kozmon, S.; Kulhánek, P.; Štepán, J.; Tvaroška, I.; Koča, J. Exploring Reaction Pathways for O-GlcNAc Transferase Catalysis. A String Method Study. J. Phys. Chem. B 2015, 119 (12), 4371−4381. (17) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. TRANSITION PATH SAMPLING: Throwing Ropes Over Rough Mountain Passes, in the Dark. Annu. Rev. Phys. Chem. 2002, 53 (1), 291−318. (18) Riplinger, C.; Neese, F. An Efficient and near Linear Scaling Pair Natural Orbital Based Local Coupled Cluster Method. J. Chem. Phys. 2013, 138 (3), 034106. (19) Liakos, D. G.; Neese, F. Is It Possible To Obtain Coupled Cluster Quality Energies at near Density Functional Theory Cost? Domain-Based Local Pair Natural Orbital Coupled Cluster vs Modern Density Functional Theory. J. Chem. Theory Comput. 2015, 11 (9), 4054−4063. (20) Fritz, T. A.; Raman, J.; Tabak, L. A. Dynamic Association between the Catalytic and Lectin Domains of Human UDP-

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00531. Figures S1−S15 (PDF)

■

Article

AUTHOR INFORMATION

Corresponding Author

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

P.J. and T.T. contributed equally. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work has been financially supported by the Ministry of Education, Youth and Sports of the Czech Republic under the project CEITEC 2020 (LQ1601) and by the Grant Agency of the Czech Republic [13-25401S]. This financial support is gratefully acknowledged. This work was supported by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project, ,IT4Innovations National Supercomputing Center − LM2015070“. Additional computational resources were provided by the CESNET LM2015042 and the CERIT Scientific Cloud LM2015085, provided under the programme “Projects of Large Research, Development, and Innovations Infrastructures”. The research has been financed by the program SASPRO (ArIDARuM, 0005/01/02 - SK) and was co-funded by the People Programme (Marie Curie Actions 7FP, grant agreement REA no. 609427 - SK) and co-financed by the Slovak Academy of Sciences. IT would like also to thank the Scientific Grant Agency of the Ministry of Education of the Slovak Republic and the Slovak Academy of Sciences (grant VEGA-02/0064/15 and 02/0024/16) for their support. 6074


Article

Journal of Chemical Theory and Computation GalNAc:Polypeptide α-N-Acetylgalactosaminyltransferase-2. J. Biol. Chem. 2006, 281 (13), 8613−8619. (21) Hagen, F. K.; Hazes, B.; Raffo, R.; deSa, D.; Tabak, L. A. Structure-Function Analysis of the UDP-N-Acetyl-D-galactosamine:PolypeptideN-Acetylgalactosaminyltransferase: ESSENTIAL RESIDUES LIE IN A PREDICTED ACTIVE SITE CLEFT RESEMBLING A LACTOSE REPRESSOR FOLD. J. Biol. Chem. 1999, 274 (10), 6797−6803. (22) Kubota, T.; Shiba, T.; Sugioka, S.; Furukawa, S.; Sawaki, H.; Kato, R.; Wakatsuki, S.; Narimatsu, H. Structural Basis of Carbohydrate Transfer Activity by Human UDP-GalNAc: Polypeptide α-N-Acetylgalactosaminyltransferase (Pp-GalNAc-T10). J. Mol. Biol. 2006, 359 (3), 708−727. (23) Stwora-Wojczyk, M. M.; Kissinger, J. C.; Spitalnik, S. L.; Wojczyk, B. S. O-Glycosylation in Toxoplasma Gondii: Identification and Analysis of a Family of UDP-GalNAc:polypeptide N-Acetylgalactosaminyltransferases. Int. J. Parasitol. 2004, 34 (3), 309−322. (24) Lee, S. S.; Hong, S. Y.; Errey, J. C.; Izumi, A.; Davies, G. J.; Davis, B. G. Mechanistic Evidence for a Front-Side, SNi-Type Reaction in a Retaining Glycosyltransferase. Nat. Chem. Biol. 2011, 7 (9), 631−638. (25) Tenno, M.; Toba, S.; Kézdy, F. J.; Elhammer, Å. P.; Kurosaka, A. Identification of Two Cysteine Residues Involved in the Binding of UDP-GalNAc to UDP-GalNAc:polypeptide N-Acetylgalactosaminyltransferase 1 (GalNAc-T1): Cys Residues in GalNAc-T1 Interact with UDP-GalNAc. Eur. J. Biochem. 2002, 269 (17), 4308−4316. (26) Tenno, M.; Saeki, A.; Elhammer, Å. P.; Kurosaka, A. Function of Conserved Aromatic Residues in the Gal/GalNAc-Glycosyltransferase Motif of UDP-GalNAc:polypeptide N-Acetylgalactosaminyltransferase 1: Structure-Function Relationship of GalNAc-T1. FEBS J. 2007, 274 (23), 6037−6045. (27) Freire, T.; FernáNdez, C.; Chalar, C.; Maizels, R. M.; Alzari, P.; Osinaga, E.; Robello, C. Characterization of a UDP-N-Acetyl-DGalactosamine:polypeptide N-Acetylgalactosaminyltransferase with an Unusual Lectin Domain from the Platyhelminth Parasite Echinococcus Granulosus. Biochem. J. 2004, 382 (2), 501−510. (28) Gómez, H.; Rojas, R.; Patel, D.; Tabak, L. A.; Lluch, J. M.; Masgrau, L. A Computational and Experimental Study of OGlycosylation. Catalysis by Human UDP-GalNAc polypeptide:GalNAc Transferase-T2. Org. Biomol. Chem. 2014, 12 (17), 2645−2655. (29) Lira-Navarrete, E.; Iglesias-Fernández, J.; Zandberg, W. F.; Compañoń , I.; Kong, Y.; Corzana, F.; Pinto, B. M.; Clausen, H.; Peregrina, J. M.; Vocadlo, D. J.; Rovira, C.; Hurtado-Guerrero, R. Substrate-Guided Front-Face Reaction Revealed by Combined Structural Snapshots and Metadynamics for the Polypeptide N -Acetylgalactosaminyltransferase 2. Angew. Chem., Int. Ed. 2014, 53 (31), 8206−8210. (30) Trnka, T.; Kozmon, S.; Tvaroška, I.; Koča, J. Stepwise Catalytic Mechanism via Short-Lived Intermediate Inferred from Combined QM/MM MERP and PES Calculations on Retaining Glycosyltransferase ppGalNAcT2. PLoS Comput. Biol. 2015, 11 (4), e1004061. (31) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Genet. 2006, 65 (3), 712−725. (32) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; GonzálezOuteiriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. GLYCAM06: A Generalizable Biomolecular Force Field. Carbohydrates. J. Comput. Chem. 2008, 29 (4), 622−655. (33) 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 (2), 926−935. (34) Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theory Comput. 2013, 9 (9), 3878−3888. (35) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B. P.; Hayik, S.; Roitberg, A. E.; Seabra, G.; Swails, J. M.;

Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Salomon-Ferrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12; University of California: San Francisco, 2012. (36) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 (8), 3684−3690. (37) CPMD; Copyright IBM Corp 1990−2008, Copyright MPI fur Festkorperforschung Stuttgart, 1997. (38) Laio, A.; VandeVondele, J.; Rothlisberger, U. A Hamiltonian Electrostatic Coupling Scheme for Hybrid Car−Parrinello Molecular Dynamics Simulations. J. Chem. Phys. 2002, 116 (16), 6941−6947. (39) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77 (18), 3865− 3868. (40) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865; Phys. Rev. Lett. 1997, 78 (7), 1396−1396. (41) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27 (15), 1787−1799. (42) Troullier, N.; Martins, J. L. Efficient Pseudopotentials for PlaneWave Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43 (3), 1993−2006. (43) Kulhánek, P.; Fuxreiter, M.; Štěpán, J.; Koča, J.; Mones, L.; Střelcová, Z.; Petřek, M. PMFLib - A Toolkit for Free Energy Calculations; Masaryk University: Brno, Czeck Republic, 2013. (44) Baerends, E. J.; Ziegler, T.; Autschbach, J.; Bashford, D.; Bérces, A.; Bickelhaupt, F. M.; Bo, C.; Boerrigter, P. M.; Cavallo, L.; Chong, D. P.; Deng, L.; Dickson, R. M.; Ellis, D. E.; van Faassen, M.; Fan, L.; Fischer, T. H.; Guerra, C. F.; Franchini, M.; Ghysels, A.; Giammona, A.; van Gisbergen, S. J. A.; Götz, A. W.; Groeneveld, J. A.; Gritsenko, O. V.; Grüning, M.; Gusarov, S.; Harris, F. E.; van den Hoek, P.; Jacob, C. R.; Jacobsen, H.; Jensen, L.; Kaminski, J. W.; van Kessel, G.; Kootstra, F.; Kovalenko, A.; Krykunov, M. V.; van Lenthe, E.; McCormack, D. A.; Michalak, A.; Mitoraj, M.; Morton, S. M.; Neugebauer, J.; Nicu, V. P.; Noodleman, L.; Osinga, V. P.; Patchkovskii, S.; Pavanello, M.; Philipsen, P. H. T.; Post, D.; Pye, C. C.; Ravenek, W.; Rodríguez, J. I.; Ros, P.; Schipper, P. R. T.; Schreckenbach, G.; Seldenthuis, J. S.; Seth, M.; Snijders, J. G.; Solà, M.; Swart, M.; Swerhone, D.; te Velde, G.; Vernooijs, P.; Versluis, L.; Visscher, L.; Visser, O.; Wang, F.; Wesolowski, T. A.; van Wezenbeek, E. M.; Wiesenekker, G.; Wolff, S. K.; Woo, T. K.; Yakovlev, A. L. ADF2013; SCM, Theoretical Chemistry, Vrije Universiteit: Amsterdam, The Netherlands, 2013. (45) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22 (9), 931−967. (46) Fonseca Guerra, C.; Snijders, J. G.; te Velde, G.; Baerends, E. J. Towards an Order-N DFT Method. Theor. Chem. Acc. 1998, 99 (6), 391−403. (47) Handy, N. C.; Cohen, A. J. Left-Right Correlation Energy. Mol. Phys. 2001, 99 (5), 403−412. (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 (15), 154104. (49) Van Lenthe, E.; Baerends, E. J. Optimized Slater-Type Basis Sets for the Elements 1−118. J. Comput. Chem. 2003, 24 (9), 1142−1156. (50) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38 (6), 3098−3100. (51) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37 (2), 785−789. 6075


Article

Journal of Chemical Theory and Computation (52) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. The Performance of a Family of Density Functional Methods. J. Chem. Phys. 1993, 98 (7), 5612−5626. (53) Russo, T. V.; Martin, R. L.; Hay, P. J. Density Functional Calculations on First-row Transition Metals. J. Chem. Phys. 1994, 101 (9), 7729−7737. (54) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98 (45), 11623−11627. (55) Perdew, J. P. Density-Functional Approximation for the Correlation Energy of the Inhomogeneous Electron Gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33 (12), 8822−8824. (56) Perdew, J. P. Erratum: Density-Functional Approximation for the Correlation Energy of the Inhomogeneous Electron Gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34 (10), 7406−7406. (57) 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 (1−3), 215−241. (58) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41 (2), 157−167. (59) Ernzerhof, M.; Scuseria, G. E. Assessment of the Perdew− Burke−Ernzerhof Exchange-Correlation Functional. J. Chem. Phys. 1999, 110 (11), 5029−5036. (60) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0Model. J. Chem. Phys. 1999, 110 (13), 6158−6170. (61) Swart, M. A New Family of Hybrid Density Functionals. Chem. Phys. Lett. 2013, 580, 166−171. (62) Tao, J.; 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 (14), 146401. (63) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Comparative Assessment of a New Nonempirical Density Functional: Molecules and Hydrogen-Bonded Complexes. J. Chem. Phys. 2003, 119 (23), 12129−12137. (64) Neese, F. The ORCA Program System. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2 (1), 73−78. (65) 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 (18), 3297−3305. (66) Weigend, F. Accurate Coulomb-Fitting Basis Sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8 (9), 1057−1065. (67) Neese, F. An Improvement of the Resolution of the Identity Approximation for the Formation of the Coulomb Matrix. J. Comput. Chem. 2003, 24 (14), 1740−1747. (68) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Efficient, Approximate and Parallel Hartree−Fock and Hybrid DFT Calculations. A “chain-of-Spheres” Algorithm for the Hartree−Fock Exchange. Chem. Phys. 2009, 356 (1−3), 98−109. (69) Karton, A.; Martin, J. M. L. Comment on: “Estimating the Hartree−Fock Limit from Finite Basis Set Calculations” [Jensen F (2005) Theor Chem Acc 113:267]. Theor. Chem. Acc. 2006, 115 (4), 330−333. (70) Truhlar, D. G. Basis-Set Extrapolation. Chem. Phys. Lett. 1998, 294 (1−3), 45−48. (71) Schwenke, D. W. The Extrapolation of One-Electron Basis Sets in Electronic Structure Calculations: How It Should Work and How It Can Be Made to Work. J. Chem. Phys. 2005, 122 (1), 014107. (72) Neese, F.; Valeev, E. F. Revisiting the Atomic Natural Orbital Approach for Basis Sets: Robust Systematic Basis Sets for Explicitly Correlated and Conventional Correlated Ab Initio Methods? J. Chem. Theory Comput. 2011, 7 (1), 33−43.

(73) Cremer, D.; Pople, J. A. General Definition of Ring Puckering Coordinates. J. Am. Chem. Soc. 1975, 97 (6), 1354−1358.

6076


MM Approaches To Elucidate ... - ACS Publications

Recommend Documents