Smooth solvation method for d-orbital semiempirical calculations of biological reactions. 1. Implementation

Altered Mechanisms for Acid-Catalyzed RNA Cleavage and Isomerization Reactions Models

RNA strand cleavage by 2'-O-transphosphorylation is catalyzed not only by numerous nucleolytic RNA enzymes (ribozymes) but also by hydroxide or hydronium ions. In experiments, both cleavage of the 5'-linked nucleoside and isomerization between 3',5'- and 2',5'-phosphodiesters occur under acidic conditions, while only the cleavage reaction is observed under basic conditions. An ab initio path-integral approach for simulating kinetic isotope effects is used to reveal the reaction mechanisms for RNA cleavage and isomerization reactions under acidic conditions. Moreover, the proposed mechanisms can also be combined through the experimental pH-rate profiles.

Abnormal Properties of Low-Dimensional Confined Water

Water molecules confined to low-dimensional spaces exhibit unusual properties compared to bulk water. For example, the alternating hydrophilic and hydrophobic nanodomains on flat silicon wafer can induce the abnormal spreading of water (contact angles near 0°) which is caused by the 2D capillary effect. Hence, exploring the physicochemical properties of confined water from the nanoscale is of great value for understanding the challenges in material science and promoting the applications of nanomaterials in the fields of mass transport, nanofluidic designing, and fuel cell. The knowledge framework of confined water can also help to better understand the complex functions of the hydration layer of biomolecules, and even trace the origin of life. In this review, the physical properties, abnormal behaviors, and functions of the confined water are mainly summarized through several common low-dimensional water formats in the fields of solid/air-water interface, nanochannel confinement, and biological hydration layer. These researches indicate that the unusual behaviors of the confined water depend strongly on the confinement size and the interaction between the molecules and confining surface. These diverse properties of confined water open a new door to materials science and may play an important role in the future development of biology.

Acceleration of Semiempirical QM/MM Methods through Message Passage Interface (MPI), Hybrid MPI/Open Multiprocessing, and Self-Consistent Field Accelerator Implementations

The strategy and implementation of scalable and efficient semiempirical (SE) QM/MM methods in CHARMM are described. The serial version of the code was first profiled to identify routines that required parallelization. Afterward, the code was parallelized and accelerated with three approaches. The first approach was the parallelization of the entire QM/MM routines, including the Fock matrix diagonalization routines, using the CHARMM message passage interface (MPI) machinery. In the second approach, two different self-consistent field (SCF) energy convergence accelerators were implemented using density and Fock matrices as targets for their extrapolations in the SCF procedure. In the third approach, the entire QM/MM and MM energy routines were accelerated by implementing the hybrid MPI/open multiprocessing (OpenMP) model in which both the task- and loop-level parallelization strategies were adopted to balance loads between different OpenMP threads. The present implementation was tested on two solvated enzyme systems (including <100 QM atoms) and an S_N2 symmetric reaction in water. The MPI version exceeded existing SE QM methods in CHARMM, which include the SCC-DFTB and SQUANTUM methods, by at least 4-fold. The use of SCF convergence accelerators further accelerated the code by ∼12-35% depending on the size of the QM region and the number of CPU cores used. Although the MPI version displayed good scalability, the performance was diminished for large numbers of MPI processes due to the overhead associated with MPI communications between nodes. This issue was partially overcome by the hybrid MPI/OpenMP approach which displayed a better scalability for a larger number of CPU cores (up to 64 CPUs in the tested systems).

VR-SCOSMO: A smooth conductor-like screening model with charge-dependent radii for modeling chemical reactions

To better represent the solvation effects observed along reaction pathways, and of ionic species in general, a charge-dependent variable-radii smooth conductor-like screening model (VR-SCOSMO) is developed. This model is implemented and parameterized with a third order density-functional tight binding quantum model, DFTB3/3OB-OPhyd, a quantum method which was developed for organic and biological compounds, utilizing a specific parameterization for phosphate hydrolysis reactions. Unlike most other applications with the DFTB3/3OB model, an auxiliary set of atomic multipoles is constructed from the underlying DFTB3 density matrix which is used to interact the solute with the solvent response surface. The resulting method is variational, produces smooth energies, and has analytic gradients. As a baseline, a conventional SCOSMO model with fixed radii is also parameterized. The SCOSMO and VR-SCOSMO models shown have comparable accuracy in reproducing neutral-molecule absolute solvation free energies; however, the VR-SCOSMO model is shown to reduce the mean unsigned errors (MUEs) of ionic compounds by half (about 2-3 kcal/mol). The VR-SCOSMO model presents similar accuracy as a charge-dependent Poisson-Boltzmann model introduced by Hou et al. [J. Chem. Theory Comput. 6, 2303 (2010)]. VR-SCOSMO is then used to examine the hydrolysis of trimethylphosphate and seven other phosphoryl transesterification reactions with different leaving groups. Two-dimensional energy landscapes are constructed for these reactions and calculated barriers are compared to those obtained from ab initio polarizable continuum calculations and experiment. Results of the VR-SCOSMO model are in good agreement in both cases, capturing the rate-limiting reaction barrier and the nature of the transition state.

Quantum, classical, and hybrid QM/MM calculations in solution: general implementation of the ddCOSMO linear scaling strategy

We present the general theory and implementation of the Conductor-like Screening Model according to the recently developed ddCOSMO paradigm. The various quantities needed to apply ddCOSMO at different levels of theory, including quantum mechanical descriptions, are discussed in detail, with a particular focus on how to compute the integrals needed to evaluate the ddCOSMO solvation energy and its derivatives. The overall computational cost of a ddCOSMO computation is then analyzed and decomposed in the various steps: the different relative weights of such contributions are then discussed for both ddCOSMO and the fastest available alternative discretization to the COSMO equations. Finally, the scaling of the cost of the various steps with respect to the size of the solute is analyzed and discussed, showing how ddCOSMO opens significantly new possibilities when cheap or hybrid molecular mechanics/quantum mechanics methods are used to describe the solute.

Ab initio path-integral calculations of kinetic and equilibrium isotope effects on base-catalyzed RNA transphosphorylation models

Detailed understandings of the reaction mechanisms of RNA catalysis in various environments can have profound importance for many applications, ranging from the design of new biotechnologies to the unraveling of the evolutionary origin of life. An integral step in the nucleolytic RNA catalysis is self-cleavage of RNA strands by 2'-O-transphosphorylation. Key to elucidating a reaction mechanism is determining the molecular structure and bonding characteristics of transition state. A direct and powerful probe of transition state is measuring isotope effects on biochemical reactions, particularly if we can reproduce isotope effect values from quantum calculations. This article significantly extends the scope of our previous joint experimental and theoretical work in examining isotope effects on enzymatic and nonenzymatic 2'-O-transphosphorylation reaction models that mimic reactions catalyzed by RNA enzymes (ribozymes), and protein enzymes such as ribonuclease A (RNase A). Native reactions are studied, as well as reactions with thio substitutions representing chemical modifications often used in experiments to probe mechanism. Here, we report and compare results from eight levels of electronic-structure calculations for constructing the potential energy surfaces in kinetic and equilibrium isotope effects (KIE and EIE) computations, including a "gold-standard" coupled-cluster level of theory [CCSD(T)]. In addition to the widely used Bigeleisen equation for estimating KIE and EIE values, internuclear anharmonicity and quantum tunneling effects were also computed using our recently developed ab initio path-integral method, that is, automated integration-free path-integral method. The results of this work establish an important set of benchmarks that serve to guide calculations of KIE and EIE for RNA catalysis.

Quantum Calculations in Solution for Large to Very Large Molecules: A New Linear Scaling QM/Continuum Approach

We present a new implementation of continuum solvation models for semiempirical Hamiltonians that allows the description of environmental effects on very large molecular systems. In this approach based on a domain decomposition strategy of the COSMO model (ddCOSMO), the solution to the COSMO equations is no longer the computational bottleneck but becomes a negligible part of the overall computation time. In this Letter, we analyze the computational impact of COSMO on the solution of the SCF equations for large to very large molecules, using semiempirical Hamiltonians, for both the new ddCOSMO implementation and the most recent, linear scaling one, based on the fast multipole method. A further analysis is on the simulation of the UV/visible spectrum of a light-harvesting pigment-protein complex. All of the results show how the new ddCOSMO algorithm paves the way to routine computations for large molecular systems in the condensed phase.

Experimental and computational analysis of the transition state for ribonuclease A-catalyzed RNA 2'-O-transphosphorylation

Enzymes function by stabilizing reaction transition states; therefore, comparison of the transition states of enzymatic and nonenzymatic model reactions can provide insight into biological catalysis. Catalysis of RNA 2'-O-transphosphorylation by ribonuclease A is proposed to involve electrostatic stabilization and acid/base catalysis, although the structure of the rate-limiting transition state is uncertain. Here, we describe coordinated kinetic isotope effect (KIE) analyses, molecular dynamics simulations, and quantum mechanical calculations to model the transition state and mechanism of RNase A. Comparison of the (18)O KIEs on the 2'O nucleophile, 5'O leaving group, and nonbridging phosphoryl oxygens for RNase A to values observed for hydronium- or hydroxide-catalyzed reactions indicate a late anionic transition state. Molecular dynamics simulations using an anionic phosphorane transition state mimic suggest that H-bonding by protonated His12 and Lys41 stabilizes the transition state by neutralizing the negative charge on the nonbridging phosphoryl oxygens. Quantum mechanical calculations consistent with the experimental KIEs indicate that expulsion of the 5'O remains an integral feature of the rate-limiting step both on and off the enzyme. Electrostatic interactions with positively charged amino acid site chains (His12/Lys41), together with proton transfer from His119, render departure of the 5'O less advanced compared with the solution reaction and stabilize charge buildup in the transition state. The ability to obtain a chemically detailed description of 2'-O-transphosphorylation transition states provides an opportunity to advance our understanding of biological catalysis significantly by determining how the catalytic modes and active site environments of phosphoryl transferases influence transition state structure.

A smooth, nonsingular, and faithful discretization scheme for polarizable continuum models: the switching/Gaussian approach

Polarizable continuum models (PCMs) are a widely used family of implicit solvent models based on reaction-field theory and boundary-element discretization of the solute/continuum interface. An often overlooked aspect of these theories is that discretization of the interface typically does not afford a continuous potential energy surface for the solute. In addition, we show that discretization can lead to numerical singularities and violations of exact variational conditions. To fix these problems, we introduce the switching/Gaussian (SWIG) method, a discretization scheme that overcomes several longstanding problems with PCMs. Our approach generalizes a procedure introduced by York and Karplus [J. Phys. Chem. A 103, 11060 (1999)], extending it beyond the conductor-like screening model. Comparison to other purportedly smooth PCM implementations reveals certain artifacts in these alternative approaches, which are avoided using the SWIG methodology. The versatility of our approach is demonstrated via geometry optimizations, vibrational frequency calculations, and molecular dynamics simulations, for solutes described using quantum mechanics and molecular mechanics.

Accurate proton affinity and gas-phase basicity values for molecules important in biocatalysis

Benchmark quantum calculations of proton affinities and gas-phase basicities of molecules relevant to biochemical processes, particularly acid/base catalysis, are presented and compared for a variety of multilevel and density functional quantum models. Included are nucleic acid bases in both keto and enol tautomeric forms, ribose in B-form and A-form sugar pucker conformations, amino acid side chains and backbone molecules, and various phosphates and phosphoranes, including thio substitutions. This work presents a high-level thermodynamic characterization of biologically relevant protonation states and provides a benchmark database for development of next-generation semiempirical and approximate density functional quantum models and parametrization of methods to predict pK(a) values and relative solvation energies.

Continuous surface charge polarizable continuum models of solvation. I. General formalism

Continuum solvation models are appealing because of the simplified yet accurate description they provide of the solvent effect on a solute, described either by quantum mechanical or classical methods. The polarizable continuum model (PCM) family of solvation models is among the most widely used, although their application has been hampered by discontinuities and singularities arising from the discretization of the integral equations at the solute-solvent interface. In this contribution we introduce a continuous surface charge (CSC) approach that leads to a smooth and robust formalism for the PCM models. We start from the scheme proposed over ten years ago by York and Karplus and we generalize it in various ways, including the extension to analytic second derivatives with respect to atomic positions. We propose an optimal discrete representation of the integral operators required for the determination of the apparent surface charge. We achieve a clear separation between "model" and "cavity" which, together with simple generalizations of modern integral codes, is all that is required for an extensible and efficient implementation of the PCM models. Following this approach we are now able to introduce solvent effects on energies, structures, and vibrational frequencies (analytical first and second derivatives with respect to atomic coordinates), magnetic properties (derivatives with respect of magnetic field using GIAOs), and in the calculation more complex properties like frequency-dependent Raman activities, vibrational circular dichroism, and Raman optical activity.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Continuous and smooth potential energy surface for conductorlike screening solvation model using fixed points with variable areas

Rigorously continuous and smooth potential energy surfaces, as well as exact analytic gradients, are obtained for a conductorlike screening solvation model (CPCM, a variant of the general COSMO) with Hartree-Fock (RHF, ROHF, UHF, and MCSCF) and density functional theory (R-DFT, RO-DFT, and U-DFT) methods using a new tessellation scheme, fixed points with variable areas (FIXPVA). In FIXPVA, spheres centered at atoms are used to define the molecular cavity and surface. The surface of each sphere is divided into 60, 240, or 960 tesserae, which have positions fixed relative to the sphere center and areas scaled by switching functions of their distances to neighboring spheres. Analytic derivatives of the positions and areas of the surface tesserae with respect to atomic coordinates can be obtained and used to evaluate the solvation energy gradients. Due to the accurate analytic gradients and smooth potential energy surface, geometry optimization processes using these methods are stable and convergent.

Extension of adaptive tree code and fast multipole methods to high angular momentum particle charge densities

The development and implementation of a tree code (TC) and fast multipole method (FMM) for the efficient, linear-scaling calculation of long-range electrostatic interactions of particle distributions with variable shape and multipole character are described. The target application of these methods are stochastic boundary molecular simulations with polarizable force fields and/or combined quantum mechanical/molecular mechanical potentials. Linear-scaling is accomplished through the adaptive decomposition of the system into a hierarchy of interacting particle sets. Two methods for effecting this decomposition are evaluated: fluc-splitting and box-splitting, for which the latter is demonstrated to be generally more accurate. In addition, a generalized termination criterion is developed that delivers optimal performance at fixed error tolerance that, in the case of quadrupole-represented Drude water, effects a speed-up by a factor of 2-3 relative to a multipole-independent termination criteria. The FMM is shown to be approximately 2-3 times faster than the TC, independent of the system size and multipole order of the particles. The TC and FMM are tested for a variety of static and polarizable water systems, and for the the 70S ribosome functional complex containing an assembly of transfer and messenger RNAs.

Electrostatic interactions in the hairpin ribozyme account for the majority of the rate acceleration without chemical participation by nucleobases

Molecular dynamics simulations using a combined quantum mechanical/molecular mechanical potential are used to determine the two-dimensional free energy profiles for the mechanism of RNA transphosphorylation in solution and catalyzed by the hairpin ribozyme. A mechanism is explored whereby the reaction proceeds without explicit chemical participation by conserved nucleobases in the active site. The ribozyme lowers the overall free energy barrier by up to 16 kcal/mol, accounting for the majority of the observed rate enhancement. The barrier reduction in this mechanism is achieved mainly by the electrostatic environment provided by the ribozyme without recruitment of active site nucleobases as acid or base catalysts. The results establish a baseline mechanism that invokes only the solvation and specific hydrogen-bonding interactions present in the ribozyme active site and provide a departure point for the exploration of alternate mechanisms where nucleobases play an active chemical role.

Charge-dependent model for many-body polarization, exchange, and dispersion interactions in hybrid quantum mechanical/molecular mechanical calculations

This work explores a new charge-dependent energy model consisting of van der Waals and polarization interactions between the quantum mechanical (QM) and molecular mechanical (MM) regions in a combined QMMM calculation. van der Waals interactions are commonly treated using empirical Lennard-Jones potentials, whose parameters are often chosen based on the QM atom type (e.g., based on hybridization or specific covalent bonding environment). This strategy for determination of QMMM nonbonding interactions becomes tedious to parametrize and lacks robust transferability. Problems occur in the study of chemical reactions where the "atom type" is a complex function of the reaction coordinate. This is particularly problematic for reactions, where atoms or localized functional groups undergo changes in charge state and hybridization. In the present work we propose a new model for nonelectrostatic nonbonded interactions in QMMM calculations that overcomes many of these problems. The model is based on a scaled overlap model for repulsive exchange and attractive dispersion interactions that is a function of atomic charge. The model is chemically significant since it properly correlates atomic size, softness, polarizability, and dispersion terms with minimal one-body parameters that are functions of the atomic charge. Tests of the model are examined for rare-gas interactions with neutral and charged atoms in order to demonstrate improved transferability. The present work provides a new framework for modeling QMMM interactions with improved accuracy and transferability.

Density functional study of the in-line mechanism of methanolysis of cyclic phosphate and thiophosphate esters in solution: insight into thio effects in RNA transesterification

Density functional calculations of thio effects on the in-line mechanism of methanolysis of ethylene phosphate (a reverse reaction model for RNA phosphate transesterification) are presented. A total of 12 reaction mechanisms are examined using the B3LYP functional with large basis sets, and the effects of solvation were treated using the PCM, CPCM, and SM5 solvation models. Single thio substitutions at all of the distinct phosphoryl oxygen positions (2', 3', 5', pro-R) and a double thio substitution at the nonbridging (pro-R/pro-S) positions were considered. Profiles for each reaction were calculated in the dianionic and monoanionic/monoprotic states, corresponding to reaction models under alkaline and nonalkaline conditions, respectively. These models provide insight into the mechanisms of RNA transesterification thio effects and serve as a set of high-level quantum data that can be used in the design of new semiempirical quantum models for hybrid quantum mechanical/molecular mechanical simulations and linear-scaling electronic structure calculations.

Improvement of semiempirical response properties with charge-dependent response density

The present work outlines a new method for treatment of charge-dependent polarizability in semiempirical quantum models for use in combined quantum-mechanical/molecular mechanical simulations of biological reactions. The method addresses a major shortcoming in the performance of conventional semiempirical models for these simulations that is tied to the use of a localized minimal atomic-orbital basis set. The present approach has the advantages that it uses a density basis that retains a set of linear-response equations, does not increase the atomic-orbital basis, and avoids the problem of artificial charge transfer and scaling of the polarizability seen in related models that allow atomic charges to fluctuate. The model introduces four new atom-based parameters and has been tested with the modified neglect of differential overlap d-orbital Hamiltonian against 1132 molecules and ions and shown to decrease the dipole moment and polarizability errors by factors of 2 and 10, respectively, with respect to density-functional results. The method performs impressively for a variety of charge states (from 2+ to 2-), and offers a potentially powerful extension in the design of next generation semiempirical quantum models for accurate simulations of highly charged biological reactions.

Smooth solvation method for d-orbital semiempirical calculations of biological reactions. 2. Application to transphosphorylation thio effects in solution

Density-functional and semiempirical quantum methods and continuum dielectric and explicit solvation models are applied to study the role of solvation on the stabilization of native and thio-substituted transphosphorylation reactions. Extensive comparison is made between results obtained from the different methods. For the semiempirical methods, explicit solvation was treated using a hybrid quantum mechanical/molecular mechanical (QM/MM) approach and the implicit solvation was treated using a recently developed smooth solvation model implemented into a d-orbital semiempirical framework (MNDO/d-SCOSMO) within CHARMM. The different quantum and solvation methods were applied to the transesterification of a 3'-ribose,5'-methyl phosphodiester that serves as a nonenzymatic model for the self-cleavage reaction catalyzed by the hammerhead and hairpin ribozymes. Thio effects were studied for a double sulfur substitution at the nonbridging phosphoryl oxygen positions. The reaction profiles of both the native and double sulfur-substituted reactions from the MNDO/d-SCOSMO calculations were similar to the QM/MM results and consistent with the experimentally observed trends. These results underscore the need for a d-orbital semiempirical representation for phosphorus and sulfur for the study of experimentally observed thio effects in enzymatic and nonenzymatic phosphoryl transfer reactions. One of the major advantages of the present approach is that it can be applied to model chemical reactions at a significantly lower computational cost than either the density-functional calculations with implicit solvation or the semiempirical QM/MM simulations with explicit solvent.

Exploring SCC-DFTB paths for mapping QM/MM reaction mechanisms

A new first-order procedure for locating transition structures (TS) that employs hybrid quantum mechanical/molecular mechanical (QM/MM) potentials has been developed. This new technique (RPATh+RESD) combines the replica path method (RPATh) and standard reaction coordinate driving (RCD) techniques in an approach that both efficiently determines reaction barriers and successfully eliminates two key weaknesses of RCD calculations (i.e., hysteresis/discontinuities in the path and the sequential nature of the RCD procedure). In addition, we have extended CHARMM's QM/MM reaction pathway methods, the RPATh and nudged elastic band (NEB) methods, to incorporate SCC-DFTB wave functions. This newly added functionality has been applied to the chorismate mutase-catalyzed interconversion of chorismate to prephenate, which is a key step in the shikimate pathway of bacteria, fungi, and other higher plants. The RPATh+RESD barrier height (DeltaE=5.7 kcal/mol) is in good agreement with previous results from full-energy surface mapping studies (Zhang, X.; Zhang, X.; Bruice, T. C. Biochemistry 2005, 44, 10443-10448). Full reaction paths were independently mapped with RPATh and NEB methods and showed good agreement with the final transition state from the RPATh+RESD "gold standard" and previous high-level QM/MM transition states (Woodcock, H. L.; Hodoscek, M.; Gilbert, T. B.; Gill, P. M. W.; Schaefer, H. F.; Brooks, B. R. J. Comput. Chem. 2007, 28, 1485-1502). The SCC-DFTB TS geometry most closely approximates the MP2/6-31+G(d) QM/MM result. However, the barrier height is underestimated and possibly points to an area for improvement in SCC-DFTB parametrization. In addition, the steepest descents (SD) minimizer for the NEB method was modified to uncouple the in-path and off-path degrees of freedom during the minimization, which significantly improved performance. The convergence behavior of the RPATh and NEB was examined for SCC-DFTB wave functions, and it was determined that, in general, both methods converge at about the same rate, although the techniques used for convergence may be different. For instance, RPATh can effectively use the adopted basis Newton-Raphson (ABNR) minimizer, where NEB seems to require a combination of SD and ABNR.

Theoretical Studies of Dissociative Phosphoryl Transfer in Interconversion of Phosphoenolpyruvate to Phosphonopyruvate: Solvent Effects, Thio Effects, and Implications for Enzymatic Reactions

The conversion of phosphoenolpyruvate (PEP) to phosphonopyruvate (P-pyr) is catalyzed by PEP mutase via a dissociative mechanism. In this work, we investigate the uncatalyzed reaction using ab initio methods, density functional theory, and the semiempirical MNDO/d method. Comparisons of geometries and relative energies of stationary points (minima and transition states) with density functional results indicate that the semiempirical method is reasonably accurate. Solvent effects are examined using implicit solvent models, including the recently extended smooth conductor-like screening model. Due to the large negative charge carried by the system, solvation is found to drastically alter the location and energy of stationary points along the dissociative reaction pathways. The influence of substituting a nonbridging phosphoryl oxygen by sulfur (thio effects) was also investigated. Implications of these results for the enzymatic reaction are discussed.