Benchmarking pK prediction methods for Lys115 in acetoacetate decarboxylase.

Three different pK prediction methods were used to calculate the pK of Lys115 in acetoacetate decarboxylase (AADase): the empirical method PROPKA, the multiconformation continuum electrostatics (MCCE) method, and the molecular dynamics/thermodynamic integration (MD/TI) method with implicit solvent. As expected, accurate pK prediction of Lys115 depends on the protonation patterns of other ionizable groups, especially the nearby Glu76. However, since the prediction methods do not explicitly sample the protonation patterns of nearby residues, this must be done manually.

Drug release by pH-responsive molecular tweezers: atomistic details from molecular modeling.

pH-responsive molecular tweezers have been proposed as an approach for targeting drug-delivery to tumors, which tend to have a lower pH than normal cells. We performed a computational study of a pH-responsive molecular tweezer using ab initio quantum chemistry in the gas-phase and molecular dynamics (MD) simulations in solution. The binding free energy in solution was calculated using steered MD.

Efficient parameterization of torsional terms for force fields.

A novel method is presented for fitting force-field dihedral angles using an ensemble of structures generated from an ab initio Monte Carlo simulation. Importance sampling is used to achieve an efficient algorithm using a low level of theory to minimize the system at each step with the dihedral angles constrained, followed by dihedral fitting using the single point energies at a higher level of theory. The resulting method is an order of magnitude more efficient than the traditional method of doing a constrained scan over each dihedral independently.

Quantum mechanics/molecular mechanics restrained electrostatic potential fitting.

We present a quantum mechanics/molecular mechanics (QM/MM) method to evaluate the partial charges of amino acid residues for use in MM potentials based on their protein environment. For each residue of interest, the nearby residues are included in the QM system while the rest of the protein is treated at the MM level of theory. After a short structural optimization, the partial charges of the central residue are fit to the electrostatic potential using the restrained electrostatic potential (RESP) method.

Efficient optimization of van der Waals parameters from bulk properties.

Due to the computational cost involved, when developing a force field for new compounds, one often avoids fitting van der Waals (vdW) terms, instead relying on a general force field based on the atom type. Here, we provide a novel approach to efficiently optimize vdW terms, based on both ab initio dimer energies and condensed phase properties. The approach avoids the computational challenges of searching the parameter space by using an extrapolation method to obtain a reliable difference quotient for the parameter derivatives based on the central difference.

Transition state analysis of enolpyruvylshikimate 3-phosphate (EPSP) synthase (AroA)-catalyzed EPSP hydrolysis.

Proton transfer to carbon atoms is a significant catalytic challenge because of the large intrinsic energetic barrier and the frequently unfavorable thermodynamics. The main catalytic challenge for enolpyruvylshikimate 3-phosphate synthase (EPSP synthase, AroA) is protonating the methylene carbon atom of phosphoenolpyruvate, or EPSP, in the reverse reaction. We performed transition state analysis using kinetic isotope effects (KIEs) on AroA-catalyzed EPSP hydrolysis, which also begins with a methylene carbon (C3) protonation, as an analog of AroA's reverse reaction.

Transition state analysis of acid-catalyzed hydrolysis of an enol ether, enolpyruvylshikimate 3-phosphate (EPSP).

Proton transfer to carbon represents a significant catalytic challenge because of the large intrinsic energetic barrier and the frequently unfavorable thermodynamics. Multiple kinetic isotope effects (KIEs) were measured for acid-catalyzed hydrolysis of the enol ether functionality of enolpyruvylshikimate 3-phosphate (EPSP) as a nonenzymatic analog of the EPSP synthase (AroA) reaction. The large solvent deuterium KIE demonstrated that protonating C3 was the rate-limiting step, and the lack of solvent hydron exchange into EPSP demonstrated that protonation was irreversible.

Automated Parametrization of AMBER Force Field Terms from Vibrational Analysis with a Focus on Functionalizing Dinuclear Zinc(II) Scaffolds.

A procedure for determining force constants that is independent of the internal redundant coordinate choice is presented. The procedure is based on solving each bond and angle term separately, using the Wilson B matrix. The method only requires a single ab initio frequency calculation at the minimum energy structure and is made available in the software "parafreq".

Computational study of the binding modes of caffeine to the adenosine A2A receptor.

Using the recently solved crystal structure of the human adenosine A(2A) receptor, we applied MM/PBSA to compare the binding modes of caffeine with those of the high-affinity selective antagonist ZM241385. MD simulations were performed in the environment of the lipid membrane bilayer. Four low-energy binding modes of caffeine-A(2A) were found, all of which had similar energies.

Interconversion study in 1,4-substituted six-membered cyclohexane-type rings. Structure and dynamics of trans-1,4-dibromo-1,4-dicyanocyclohexane.

Cyclohexane is an extremely flexible molecule that oscillates, at room temperature, between two clearly distinct and extreme conformations that cannot be distinguished at room temperature; so much so that the NMR spectrum is a single line that includes all 12 protons be they axial or equatorial. This raises the interesting question as to what happens when there are equal substituents at the 1 and 4 carbon atoms of the ring. Therefore substitution in the 1,4-positions in the cyclohexane ring has been the subject of considerable interest because some form of interconversion between extreme conformations could lead to the existence of a rather unusual behavior.

A parameterized, continuum electrostatic model for predicting protein pKa values.

Recognizing the limits of trying to achieve chemical accuracy for pK(a) calculations with a purely electrostatic model, we include empirical corrections into the Poisson-Boltzmann solver macroscopic electrostatics with atomic detail (Bashford, Biochemistry 1990;29:10219-10225), to improve the reliability and accuracy of the model. The total number of parameters is kept to a minimum to maximize the robustness of the model for compounds outside of the fitting dataset. The parameters are based on: (a) the electrostatic interaction between functional groups close to the titratable site, (b) the electrostatic work required to desolvate the residue, and (c) the site-to-site interactions.

Empirical prediction of protein pKa values with residue mutation.

A fast, empirical method, Mut-pKa, is presented for predicting the pKa values of ionizable residues in proteins based on mutation. The method compares the effect of mutating each residue that may act as a hydrogen bond donor or acceptor for the ionizable residue. The energetic effect of each type of mutation, along with a desolvation measure and the overall background charge, is fit against pKa data for histidine and carboxyl residues.

Practical calculation of molecular acidity with the aid of a reference molecule.

A set of linear free energy models are presented for determining the pK(a) values of amines, alcohols, and carboxylic acids. Models are determined from a series of pK(a) predictors, taken both from traditional natural atomic orbital analysis (NAO) and from a novel approach introduced here of using a reference molecule: an ammonium ion for amines and a hydrogen sulfide molecule for alcohols and carboxylic acids. Using these reference molecules, we calculate the barrier to proton transfer and show that a number of properties associated with the transition state are correlated with the pK(a).

Quantum mechanics/molecular mechanics strategies for docking pose refinement: distinguishing between binders and decoys in cytochrome C peroxidase.

We investigate the effect of systematically applying molecular dynamics (MD) and quantum mechanics/molecular mechanics (QM/MM) to docked poses in an attempt to improve the correspondence between theoretical prediction and experimental observation. The proposed scheme involves running a short time scale MD simulation on a docked ligand pose (and any known structurally important crystal structure waters in the active site), followed by QM/MM minimization. Both of these steps are relatively fast for moderately sized ligands; longer time scale MD involving the protein is not found to improve the results.

Quasi-Newton parallel geometry optimization methods.

Algorithms for parallel unconstrained minimization of molecular systems are examined. The overall framework of minimization is the same except for the choice of directions for updating the quasi-Newton Hessian. Ideally these directions are chosen so the updated Hessian gives steps that are same as using the Newton method.

Methods for finding transition states on reduced potential energy surfaces.

Three new algorithms are presented for determining transition state (TS) structures on the reduced potential energy surface, that is, for problems in which a few important degrees of freedom can be isolated. All three methods use constrained optimization to rapidly find the TS without an initial Hessian evaluation. The algorithms highlight how efficiently the TS can be located on a reduced surface, where the rest of the degrees of freedom are minimized.

Dual Grid Methods for Finding the Reaction Path on Reduced Potential Energy Surfaces.

Two new algorithms are presented for determining the minimum energy reaction path (MEP) on the reduced potential energy surface (RPES) starting with only the reactant. These approaches are based on concepts from the fast marching method (FMM), which expands points outward as a wavefront on a multidimensional grid from the reactant until the product is reached. The MEP is then traced backward to the reactant.

Moving least-squares enhanced Shepard interpolation for the fast marching and string methods.

The number of the potential energy calculations required by the quadratic string method (QSM), and the fast marching method (FMM) is significantly reduced by using Shepard interpolation, with a moving least squares to fit the higher-order derivatives of the potential. The derivatives of the potential are fitted up to fifth order. With an error estimate for the interpolated values, this moving least squares enhanced Shepard interpolation scheme drastically reduces the number of potential energy calculations in FMM, often by up 80%.

Linear-scaling quantum calculations using non-orthogonal localized molecular orbitals.

An absolute energy minimum variational principle is used for carrying out linear scaling calculations with non-orthogonal localized orbitals. Compared with results based on orthogonal localized molecular orbitals, the method is shown to give significantly more accurate results when the localized molecular orbitals are allowed to be non-orthogonal. This is made possible by introducing a second minimization for approximating the inverse overlap matrix.

Quantum mechanics/molecular mechanics minimum free-energy path for accurate reaction energetics in solution and enzymes: sequential sampling and optimization on the potential of mean force surface.

To accurately determine the reaction path and its energetics for enzymatic and solution-phase reactions, we present a sequential sampling and optimization approach that greatly enhances the efficiency of the ab initio quantum mechanics/molecular mechanics minimum free-energy path (QM/MM-MFEP) method. In the QM/MM-MFEP method, the thermodynamics of a complex reaction system is described by the potential of mean force (PMF) surface of the quantum mechanical (QM) subsystem with a small number of degrees of freedom, somewhat like describing a reaction process in the gas phase. The main computational cost of the QM/MM-MFEP method comes from the statistical sampling of conformations of the molecular mechanical (MM) subsystem required for the calculation of the QM PMF and its gradient.

Sequential quadratic programming method for determining the minimum energy path.

A new method, referred to as the sequential quadratic programming method, is presented for determining minimum energy paths. The method is based on minimizing the points representing the path in the subspace perpendicular to the tangent of the path while using a penalty term to prevent kinks from forming. Rather than taking one full step, the minimization is divided into a number of sequential steps on an approximate quadratic surface.

Automatic integration of the reaction path using diagonally implicit Runge-Kutta methods.

The diagonally implicit Runge-Kutta framework is shown to be a general form for constructing stable, efficient steepest descent reaction path integrators, of any order. With this framework tolerance driven, adaptive step-size methods can be constructed by embedding methods to obtain error estimates of each step without additional computational cost. There are many embedded and nonembedded, diagonally implicit Runge-Kutta methods available from the numerical analysis literature and these are reviewed for orders two, three, and four.

A combined explicit-implicit method for high accuracy reaction path integration.

We present the use of an optimal combined explicit-implicit method for following the reaction path to high accuracy. This is in contrast to most purely implicit reaction path integration algorithms, which are only efficient on stiff ordinary differential equations. The defining equation for the reaction path is considered to be stiff, however, we show here that the reaction path is not uniformly stiff and instead is only stiff near stationary points.

Quadratic string method for determining the minimum-energy path based on multiobjective optimization.

Based on a multiobjective optimization framework, we develop a new quadratic string method for finding the minimum-energy path. In the method, each point on the minimum-energy path is minimized by integration in the descent direction perpendicular to path. Each local integration is done on a quadratic surface approximated by a damped Broyden-Fletcher-Goldfarb-Shanno updated Hessian, allowing the algorithm to take many steps between energy and gradient calls.