Prediction of pKa shifts in proteins using a combination of molecular mechanical and continuum solvent calculations

The prediction of pKa shifts of ionizable groups in proteins is of great relevance for a number of important biological phenomena. We present an implementation of the MM-GBSA approach, which combines molecular mechanical (MM) and generalized Born (GB) continuum solvent energy terms, to the calculation of pKa values of a panel of nine proteins, including 69 individual comparisons with experiment. While applied so far mainly to the calculation of biomolecular binding free energies, we show that this method can also be used for the estimation of protein pKa shifts, with an accuracy around 1 pKa unit, even for strongly shifted residues. Our analysis reveals that the nonelectrostatic terms that are part of the MM-GBSA free energy expression are important contributors to improved prediction accuracy. This suggests that most of the previous approaches that focus only on electrostatic interactions could be improved by adding other nonpolar energy terms to their free energy expression. Interestingly, our method yields best accuracy at protein dielectric constants of epsilonint = 2-4, which is in contrast to previous approaches that peak at higher epsilonint > or = 8. An important component of our procedure is an intermediate minimization step of each protonation state involving different rotamers and tautomers as a way to explicitly model protein relaxation upon (de)protonation.     

Theoretical pKa calculations of proteins; the tyrosine and lysine residues of b-elicitin

Elicitins are small proteins that are secreted by plant pathogenic fungi. In this work we have used a computer program that utilizes the boundary element method for heterogeneous dielectrics with ionic strength to calculate the pK _a of all titrating groups in the 98-residue protein β-cryptogein. Our results are in reasonable agreement with the experimentally determined pK _a values for the Tyr residues in the protein. We find that the functionally important Lys13 residue has a normal pK _a of 10.3. Our work also shows that there is no direct correlation between the exposure of an amino acid sidechain and its pK _a. Received: 24 April 1998 / Accepted: 4 August 1998 / Published online: 11 November 1998     

Accurate pKa determination for a heterogeneous group of organic molecules.

Single-molecule studies that allow to compute pKa values, proton affinities (gas-phase acidity/basicity) and the electrostatic energy of solvation have been performed for a heterogeneous set of 26 organic compounds. Quantum mechanical density functional theory (DFT) using the Becke-half&half and B3LYP functionals on optimized molecular geometries have been carried out to investigate the energetics of gas-phase protonation. The electrostatic contribution to the solvation energies of protonated and deprotonated compounds were calculated by solving the Poisson equation using atomic charges generated by fitting the electrostatic potential derived from the molecular wave functions in vacuum. The combination of gas-phase and electrostatic solvation energies by means of the thermodynamic cycle enabled us to compute pKa values for the 26 compounds, which cover six distinct chemical groups (carboxylic acids, benzoic acids, phenols, imides, pyridines and imidazoles). The computational procedure for determining pKa values is accurate and transferable with a root-mean-square deviation of 0.53 and 0.57 pKa units and a maximum error of 1.0 pKa and 1.3 pKa units for Becke-half&half and B3LYP DFT functionals, respectively.     

Evidence that there is a future for semiempirical molecular orbital calculations

The frequency of use and discussion of semiempirical and ab initio software is traced with bibliometric data from the Current Journals of the American Chemistry Society (JCACS) database, which has complete papers published in 19 journals in various fields of chemistry. Not only is the use of all types of computational chemistry software increasing apace, but also the number of papers mentioning semiempirical programs as a percentage of all papers mentioning quantum chemistry programs is growing. This percentage referring to semiempirical molecular orbital software grew from 29% in 1989 to 34% in 1993.     

Determination of pKa values for (E)-2-hydroxy-5-(aryldiazenyl) benzaldehydes in dimethyl sulfoxide: Cyclic voltammetry and density functional theory calculations

The (E)-2-hydroxy-5-(aryldiazenyl) benzaldehydes (azo dyes 1–4 ) were synthesized in high purity. As they are insoluble in water, the usual analytical methods cannot be utilized to determine their pK_a values. Cyclic voltammetry was experimentally used to determine their pKa values in DMSO solvent. In addition, computational methods and a conductor-like screening model (COSMO) were used to calculate the solvent effect. , , K_exchange, and pK_a values were estimated for the azo dyes being studied using the BP86, TPSS, B3LYP, PBE0, TPSSh, and PW6B95 density functionals in def2-TZVP basis sets. The obtained mean absolute deviations (MADs) indicate that the results of BP86, PBE0, and PW6B95 functionals are in good agreement with experimental values.     

Study on the behavior of energy convergence in ab initio crystal orbital calculations

Summary This article studies the dependence on the cutoff scheme of ab initio crystal orbital calculations with no long-range correction. We have thoroughly studied the Namur cutoff and cell-wise cutoff schemes through calculations of polyethylene and LiH chains. The Namur cutoff gives the fastest energy convergence with respect to the number of neighbors (N ₀). The energy convergence behavior with respect to N ₀ depends on the basis set. The Namur cutoff shows the fastest convergence with the STO-3G basis set, intermediate convergence with the MINI basis set, and the slowest convergence with the (7s4p/3s) basis set. The cell-wise cutoff shows exactly the reverse order of the Namur cutoff. The Namur cutoff destroys the translational symmetry. Both the Namur cutoff and cell-wise cutoff schemes introduce slight asymmetry on the two equivalent C-C bonds of polyethylene when calculating with a C₂H₄ unit cell. The asymmetry with the Namur cutoff can be made to disappear by increasing N ₀ a little. The calculations on two different unit-cell structures of trans-polyacetylene show the effect of the cutoff scheme on the total energy. Only the symmetric cutoff energies are the same. Disagreement related to the Namur cutoff disappears at N ₀ = 20, however, that related to the cell-wise and modified symmetric cutoff schemes remains at N ₀ 20. The optimized geometry and vibrational frequency are not as sensitive to the cutoff method except with the symmetric cutoff. A compilation of all results shows that the Namur cutoff is the superior cutoff scheme when calculating the insulator using the minimal basis set, especially the STO-3G basis set.     

Using polarizable POSSIM force field and fuzzy‐border continuum solvent model to calculate pKa shifts of protein residues

Our Fuzzy‐Border (FB) continuum solvent model has been extended and modified to produce hydration parameters for small molecules using POlarizable Simulations Second‐order Interaction Model (POSSIM) framework with an average error of 0.136 kcal/mol. It was then used to compute pK _a shifts for carboxylic and basic residues of the turkey ovomucoid third domain (OMTKY3) protein. The average unsigned errors in the acid and base pK _a values were 0.37 and 0.4 pH units, respectively, versus 0.58 and 0.7 pH units as calculated with a previous version of polarizable protein force field and Poisson Boltzmann continuum solvent. This POSSIM/FB result is produced with explicit refitting of the hydration parameters to the pK _a values of the carboxylic and basic residues of the OMTKY3 protein; thus, the values of the acidity constants can be viewed as additional fitting target data. In addition to calculating pK _a shifts for the OMTKY3 residues, we have studied aspartic acid residues of Rnase Sa. This was done without any further refitting of the parameters and agreement with the experimental pK _a values is within an average unsigned error of 0.65 pH units. This result included the Asp79 residue that is buried and thus has a high experimental pK _a value of 7.37 units. Thus, the presented model is capable or reproducing pK _a results for residues in an environment that is significantly different from the solvated protein surface used in the fitting. Therefore, the POSSIM force field and the FB continuum solvent parameters have been demonstrated to be sufficiently robust and transferable. © 2016 Wiley Periodicals, Inc.     

Use of semiempirical molecular orbital calculations for the evaluation of force fields

Various methods, employing molecular orbital calculations of varying approximations, for evaluation of force fields of polyatomic molecules have been reviewed. Applications ofcndo/force method for the force field calculations are specially dealt with in detail because of its ease of operation and being economically more viable in terms of computer time. The calculated C=O stretching force constants for a series of organic molecules are shown to have linear relationship with substituent constants.     

Theoretical calculation of the pKa values of some drugs in aqueous solution

In this work, calculations of pK_a values have been performed on benzoic acid and its para‐substituted derivatives and some drugs by using Gaussian 98 software package. Gas‐phase energies were calculated with HF/6‐31 G** and B3LYP/6‐31 G** levels of theory. Free energies of solvation have been computed using the polarizable continuum model (PCM), conductor‐like PCM (CPCM), and the integral equation formalism‐PCM at the same levels which have been used for geometry determination in the gas‐phase. The results that show the calculated pK_a values using the B3LYP are better than those using the corresponding HF. In comparison to the other models, the results obtained indicate that the PCM model is a suitable solvation model for calculating pK_a values. For the investigated compounds, a good agreement between the experimental and the calculated pK_a values was also observed. © 2011 Wiley Periodicals, Inc. Int J Quantum Chem, 2011     

Sparse matrix multiplications for linear scaling electronic structure calculations in an atom-centered basis set using multiatom blocks

A sparse matrix multiplication scheme with multiatom blocks is reported, a tool that can be very useful for developing linear-scaling methods with atom-centered basis functions. Compared to conventional element-by-element sparse matrix multiplication schemes, efficiency is gained by the use of the highly optimized basic linear algebra subroutines (BLAS). However, some sparsity is lost in the multiatom blocking scheme because these matrix blocks will in general contain negligible elements. As a result, an optimal block size that minimizes the CPU time by balancing these two effects is recovered. In calculations on linear alkanes, polyglycines, estane polymers, and water clusters the optimal block size is found to be between 40 and 100 basis functions, where about 55-75% of the machine peak performance was achieved on an IBM RS6000 workstation. In these calculations, the blocked sparse matrix multiplications can be 10 times faster than a standard element-by-element sparse matrix package.     

Calculating pKa values for substituted phenols and hydration energies for other compounds with the first‐order fuzzy‐border continuum solvation model

We have computed pK_a values for 11 substituted phenol compounds using the continuum Fuzzy‐Border (FB) solvation model. Hydration energies for 40 other compounds, including alkanes, alkenes, alkynes, ketones, amines, alcohols, ethers, aromatics, amides, heterocycles, thiols, sulfides, and acids have been calculated. The overall average unsigned error in the calculated acidity constant values was equal to 0.41 pH units and the average error in the solvation energies was 0.076 kcal/mol. We have also reproduced pK_a values of propanoic and butanoic acids within about 0.1 pH units from the experimental values by fitting the solvation parameters for carboxylate ion carbon and oxygen atoms. The FB model combines two distinguishing features. First, it limits the amount of noise which is common in numerical treatment of continuum solvation models by using fixed‐position grid points. Second, it uses either second‐ or first‐order approximation for the solvent polarization, depending on a particular implementation. These approximations are similar to those used for solute and explicit solvent fast polarization treatment which we developed previously. This article describes results of using the first‐order technique. This approximation places the presented methodology between the Generalized Born and Poisson‐Boltzmann continuum solvation models with respect to their accuracy of reproducing the many‐body effects in modeling a continuum solvent. © 2012 Wiley Periodicals, Inc.     

Peptide models XXXI. Conformational properties of hydrophobic residues shaping the core of proteins. An ab initio study of N‐formyl‐L‐valinamide and N‐formyl‐L‐phenylalaninamide

Employing introductory (3‐21G RHF) and medium‐size (6‐311++G** B3LYP) ab initio calculations, complete conformational libraries, containing as many as 27 conformers, have been determined for diamide model systems incorporating the amino acids valine (Val) and phenylalanine (Phe). Conformational and energetic properties of these libraries were analyzed. For example, significant correlation was found between relative energies from 6‐311++G** B3LYP and single‐point B3LYP/6‐311++G**//RHF/3‐21G calculations. Comparison of populations of molecular conformations of hydrophobic aromatic and nonaromatic residues, based on their ab initiorelative energies, with their natural abundance indicates that, at least for the hydrophobic core of proteins, the conformations of Val (Ile, Leu) and Phe (Tyr, Trp) are controlled by the local energetic preferences of the respective amino acids. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 732–751, 2001     

Colle-Salvetti-type correction for electron-nucleus correlation in the nuclear orbital plus molecular orbital theory

A Colle-Salvetti (CS)-type electron-nucleus correction in the nuclear orbital plus molecular orbital theory is proposed. The CS-type correction is designed to satisfy the cusp condition for the electron-nucleus interaction. Since the CS-type correction is expressed in terms of the electron and nucleus densities, its evaluation is computationally feasible. Numerical assessment confirms that the CS-type correction performs well for the small G2 set.     

Ab initio calculations of pKa values of some organic acids in aqueous solution

The acidity of different classes of organic compounds in aqueous solution has been calculated. The calculations are carried out at the SCF level with inclusion of entropic and thermochemical correction to yield free energies of dissociations.

The polarized continuum model is used to describe the solvent. The model furnishes pK_a values in relatively good agreement with experimental data. Scaling different parts of solvation energies provides a significant improvement in results and signifies the importance of balance of individual contributions from electrostatic, cavity, dispersion and repulsion interactions.     

Kinetics and molecular orbital calculations of desulfonation of strong acid cation-exchange resin

This paper deals with the desulfonation properties of some strong acid cation-exchange resins. The sulfate concentration in solution is continuously increased when a strong acid cation-exchange resin is mixed with water. The leaching of sulfate results from the desulfonation of the fixed group, and the amount of leached sulfate depends on the counter ion charge, the crosslinking degree and the exchanger matrix. The effects of the counter ion charge on the desulfonation rate suggested that the counter ion induces the nucleophilic attack of a water molecule on the sulfo group. This interpretation was supported by semiempirical molecular orbital calculations for the C₆H₅SO₃⁻M^m+ (M^m+ = Na⁺, Mg²⁺ and Al³⁺) systems, and the transition state of the Na⁺ system was successfully predicted by DFT calculations. The crosslinking degree influenced the desulfonation rate, which can be related to the decreasing hydration number of each counter ion in the resin phase with the increasing crosslinking degree. Furthermore, different exchanger matrices produced the differences in the rates, which may be derived from the electron-density donation from the exchanger matrix to the sulfo group. The desulfonation is governed by the electron-density of the sulfur atom and the water activity in the solid phase.     

A molecular orbital rationalization of ligand effects in N2 activation

Molecular orbital theory has been used to study a series of [(micro-N2){ML3}2] complexes as models for dinitrogen activation, with M=Mo, Ta, W, Re and L=NH2, PH2, AsH2, SbH2 and N(BH2)2. The main aims of this study have been to provide a thorough electronic analysis of the complexes and to extend previous work involving molecular orbital analyses. Molecular orbital diagrams have been used to rationalize why for L=NH2 ligand rotation is important for the singlet state but not the triplet, to confirm the effect of ligand pi donation, and to rationalize the importance of the metal d-electron configuration. The outcomes of this study will assist with a more in-depth understanding of the electronic basis for N2 activation and allow clearer predictions to be made about the structure and multiplicity of systems involved in transition-metal catalysis.     

Semi-empirical (PM3 and AM1) and ab initio molecular orbital calculations of 1,2,4-oxadiazoles, 4,5-dihydro-1,2,4-oxadiazoles and 4,4-di-n-butyl-2-phenylbenzo-1,3-oxazine

Molecular orbital calculations using semi-empirical (PM3 and AM1) and ab initio (HF/6-31G) types have been carried out on several 3,5-disubstituted 1,2,4-oxadiazoles 1a–d, 5-n-butyl-3,5-diaryl-4,5-dihydro-1,2,4-oxadiazoles 2a–d, and 4,4-di-n-butyl-2-phenylbenzo-1,3-oxazine 3. A comparison of the results by the two computational procedures has been made. Transformation of the oxadiazole ring to 4,5-dihydro-1,2,4-oxadiazole having both aryl and n-butyl groups at C-5 exhibited interesting conformational features. Also, examination of 1,3-oxazine 3 gave an idea about the structure of this compound. The rotational barrier of each phenyl group in 1a and 1d has been calculated using the ab initio method HF/6-31G(d).     

Using the locality of the small-component density in molecular Dirac–Hartree-Fock calculations

A novel scheme to improve the computational efficiency of the Dirac–Hartree-Fock method was implemented and tested in different model systems. The method uses a one-center approximation to remove all multicenter electron-repulsion integrals over the small-component basis. In all cases we found the associated errors to be below chemical accuracy, which makes the method suitable for routine application to molecules that contain heavy elements. Received: 21 December 2001 / Accepted: 11 February 2002 / Published online: 3 May 2002     

Generalized hybrid orbital for the treatment of boundary atoms in combined quantum mechanical and molecular mechanical calculations using the semiempirical parameterized model 3 method

The application of combined quantum mechanical (QM) and molecular mechanical methods to large molecular systems requires an adequate treatment of the boundary between the two approaches. In this article, we extend the generalized hybrid orbital (GHO) method to the semiempirical parameterized model 3 (PM3) Hamiltonian combined with the CHARMM force field. The GHO method makes use of four hybrid orbitals, one of which is included in the QM region in self-consistent field optimization and three are treated as auxiliary orbitals that do not participate in the QM optimization, but they provide an effective electric field for interactions. An important feature of the GHO method is that the semiempirical parameters for the boundary atom are transferable, and these parameters have been developed for a carbon boundary atom consistent with the PM3 model. The combined GHO-PM3/CHARMM model has been tested on molecular geometry and proton affinity for a series of organic compounds.Acknowledgement We thank the National Institutes of Health for support of this research.Contribution to the Jacopo Tomasi Honorary Issue