A Combined QM/MM Poisson−Boltzmann Approach

Generalized energy-based fragmentation approach for calculations of solvation energies of large systems

A generalized energy-based fragmentation (GEBF) approach has been combined with a universal solvation model based on solute electron density (SMD) to compute the solvation energies of general large systems (such as protein molecules) in solutions. In the GEBF-SMD method, the solvation energy of a target system could be combined by the corresponding solvation energies of various subsystems, each of which is embedded in the background point charges and surface charges on the surface of solute cavity at the positions of its atoms and neighbouring atoms outside of the subsystem. Our results show that the GEBF-SMD model could reproduce the conventional SMD solvation energies quite well for various proteins in solutions, and could significantly reduce the computational costs for the SMD calculations of large proteins. In addition, the GEBF-SMD approach is almost independent of the basis sets and the types of solvents (including protic, polar, and nonpolar ones). Also, the GEBF-SMD approach could reproduce the relative energies of various conformers of large systems in solutions. Therefore, the GEBF-SMD method is expected to be applicable for computing the solvation energies of a broad range of large systems.

Importance of Lattice Constants in QM/MM Calculations on Metal-Organic Frameworks

Metal-organic frameworks (MOFs) are crystalline materials with novel physical and chemical properties. Computational simulations have become powerful complements to experiment for understanding catalysis in MOFs and developing new MOFs and their applications. However, due to their relatively large and complex structures, MOFs can be burdensome for fully quantum mechanical calculations. A combined quantum mechanical and molecular mechanical (QM/MM) method that combines the accuracy of fully quantum mechanical methods and the efficiency of MM methods is therefore attractive. In this study, we employ a QM/MM method for the study of two classes of chemical process in a MOF: the conversion of reaction intermediates in an Ir-containing borylation catalyst supported on MOF UiO-67 and the diffusion of a diborylated methane molecule in the pristine UiO-67 framework. We compare the QM/MM results with full-quantum mechanical results on large systems to validate the accuracy of the applied QM/MM method. In the first case, we consider a model of the entire system by partitioning it into subsystems that interact covalently, and in the second case the subsystem interaction is mainly steric. We observe that the QM/MM results agree with the full-quantum mechanical results within an average of 4 kcal/mol in the first case with strong electronic interactions and within an average of 3 kcal/mol in the case with only noncovalent interactions. An important lesson learned from the present study is that the quantitative results are very sensitive to the lattice constants predicted by the MM method used in the QM/MM calculations.

Incorporating QM and solvation into docking for applications to GPCR targets

A great number of GPCR crystal structures have been solved in recent years, enabling GPCR-targeted drug discovery using structure-based approaches such as docking. GPCRs generally have wide and open entrances to the binding sites, which render the binding sites readily accessible to solvent. GPCRs are also populated with hydrophilic residues in the extracellular regions. Thus, including solvent and polarization effects can be important for accurate GPCR docking. To test this hypothesis, a new docking protocol which incorporates quantum mechanical/molecular mechanical (QM/MM) calculations along with an implicit solvent model is developed. The new docking method treats the ligands and the protein residues in the binding sites as QM regions and performs QM/MM calculations with implicit solvent. The results of a test on all solved GPCR cocrystals show a significant improvement over the conventional docking method.

A General Boundary Potential for Hybrid QM/MM Simulations of Solvated Biomolecular Systems

We present a general boundary potential for the efficient and accurate evaluation of electrostatic interactions in hybrid quantum mechanical/molecular mechanical (QM/MM) approaches called solvated macromolecule boundary potential (SMBP), which is designed for QM/MM calculations with any kind of QM method. The SMBP targets QM/MM single-point energy calculations and geometry optimizations. In the SMBP scheme, the outer solvent and macromolecule region is described by a boundary potential obtained with the use of Poisson-Boltzmann calculations (treating the bulk solvent as a dielectric continuum). In the QM calculations, the SMBP is represented by virtual point charges on a surface enclosing the explicitly treated inner region. These charges and their interactions with the QM density are determined through a self-consistent reaction field procedure. The accuracy of the SMBP is evaluated on three diverse test systems: the intramolecular proton transfer of glycine in water, the hydroxylation reaction in p-hydroxybenzoate hydroxylase, and the spin state energy splittings in the pentacoordinated ferric complex of cytochrome P450cam. In the case of solvated glycine, application of the SMBP turns out to be problematic since analogous QM/MM/SMBP and full QM/MM geometry optimizations lead to different close-lying local minima. In both enzymes, the SMBP performs very well and closely reproduces the results from full QM/MM optimizations of these more rigid test systems. Starting from optimized QM/MM/SMBP structures along a reaction path, one can apply the previously implemented generalized solvent boundary potential (GSBP) to sample over MM phase space in QM/MM free energy calculations within the framework of free energy perturbation theory. This reduces the overall computational costs of sampling by 1 order of magnitude while maintaining good accuracy. The combined use of SMBP and GSBP thus allows for efficient QM/MM free energy studies of enzymes.

A "Stepping Stone" Approach for Obtaining Quantum Free Energies of Hydration

We present a method which uses DFT (quantum, QM) calculations to improve free energies of binding computed with classical force fields (classical, MM). To overcome the incomplete overlap of configurational spaces between MM and QM, we use a hybrid Monte Carlo approach to generate quickly correct ensembles of structures of intermediate states between a MM and a QM/MM description, hence taking into account a great fraction of the electronic polarization of the quantum system, while being able to use thermodynamic integration to compute the free energy of transition between the MM and QM/MM. Then, we perform a final transition from QM/MM to full QM using a one-step free energy perturbation approach. By using QM/MM as a stepping stone toward the full QM description, we find very small convergence errors (<1 kJ/mol) in the transition to full QM. We apply this method to compute hydration free energies, and we obtain consistent improvements over the MM values for all molecules we used in this study. This approach requires large-scale DFT calculations as the full QM systems involved the ligands and all waters in their simulation cells, so the linear-scaling DFT code ONETEP was used for these calculations.

Aqueous solvation of polyalanine α-helices with specific water molecules and with the CPCM and SM5.2 aqueous continuum models using density functional theory

We present density functional theory (DFT) calculations at the X3LYP/D95(d,p) level on the solvation of polyalanine α-helices in water. The study includes the effects of discrete water molecules and the CPCM and AMSOL SM5.2 solvent continuum model both separately and in combination. We find that individual water molecules cooperatively hydrogen-bond to both the C- and N-termini of the helix, which results in increases in the dipole moment of the helix/water complex to more than the vector sum of their individual dipole moments. These waters are found to be more stable than in bulk solvent. On the other hand, individual water molecules that interact with the backbone lower the dipole moment of the helix/water complex to below that of the helix itself. Small clusters of waters at the termini increase the dipole moments of the helix/water aggregates, but the effect diminishes as more waters are added. We discuss the somewhat complex behavior of the helix with the discrete waters in the continuum models.

An extrapolation method for computing protein solvation energies based on density fragmentation of a graphical surface tessellation

Modeling chemical events inside proteins often require the incorporation of solvent effects via continuum polarizable models. One of these approaches is based on the assumption that the interface between solute and solvent acts as a conductor. Image charges are added on the molecular surface to satisfy the appropriate conductor boundary conditions in the presence of solute charges. As in the case of other polarizable continuum models that are based on surface tessellation, the simplest implementation of this approach is often limited to several hundred atoms due to a matrix inversion, which scales as the cube of the number or tesserae. For larger systems, approaches that use iterative matrix solvers coupled to fast summation methods must be used. In the present work, we develop a self-consistent approach to obtain conductor-like screening charges suitable for applications in proteins. The approach is based on a density fragmentation of a graphical surface tessellation. This method, although approximate, provides a straightforward scheme of parallelization, which can in principle be added to existing linear scaling implementations of conductor-like models. We implement this method in conjunction with a fixed charge model for the protein, as well as with a moving domain QM/MM description of the protein. In the latter case, the overall result leads to a charge distribution within the protein determined by self-polarization and polarization due to solvent.

A Mixed QM/MM Scoring Function to Predict Protein-Ligand Binding Affinity

Computational methods for predicting protein-ligand binding free energy continue to be popular as a potential cost-cutting method in the drug discovery process. However, accurate predictions are often difficult to make as estimates must be made for certain electronic and entropic terms in conventional force field based scoring functions. Mixed quantum mechanics/molecular mechanics (QM/MM) methods allow electronic effects for a small region of the protein to be calculated, treating the remaining atoms as a fixed charge background for the active site. Such a semi-empirical QM/MM scoring function has been implemented in AMBER using DivCon and tested on a set of 23 metalloprotein-ligand complexes, where QM/MM methods provide a particular advantage in the modeling of the metal ion. The binding affinity of this set of proteins can be calculated with an R(2) of 0.64 and a standard deviation of 1.88 kcal/mol without fitting and 0.71 and a standard deviation of 1.69 kcal/mol with fitted weighting of the individual scoring terms. In this study we explore using various methods to calculate terms in the binding free energy equation, including entropy estimates and minimization standards. From these studies we found that using the rotational bond estimate to ligand entropy results in a reasonable R(2) of 0.63 without fitting. We also found that using the ESCF energy of the proteins without minimization resulted in an R(2) of 0.57, when using the rotatable bond entropy estimate.

On the applicability of fragmentation methods to conjugated pi systems within density functional framework

For the accurate ab initio treatment of large molecular systems, linear scaling methods (LSMs) have been devised and successfully applied to covalently bonded systems as well as to those involving weak intra/intermolecular bonds. Very few attempts to apply LSM to highly conjugated molecules, especially to two-dimensional systems, have so far been reported in the literature. The present article examines the applicability of a LSM, viz., molecular tailoring approach (MTA), to pi-conjugated systems within density functional theory. A few test cases within second order Møller-Plesset framework are also reported. MTA is applied to some one-dimensional pi-conjugated molecules, for which the difference between MTA energy and actual energy is found out to be less than 1 mhartree and also reduced computation time as well as hardware requirements. The method is also extended to some small/medium-sized two-dimensional pi-conjugated molecules by developing a systematic algorithm for tailoring such systems. However, for such systems, although the energies are in error by a few millihartrees, gradients are found to match reasonably well their actual counterparts. Hence, geometry optimization of these systems within MTA framework is attempted. The geometries thus generated are found to be in good agreement with their actual counterparts, with the actual single point energies matching within 1 mhartree, along with reduced computational effort. These results point toward the potential applicability of MTA to large two- and three-dimensional pi-conjugated systems.

Polarization and polarizability assessed by protein amide acidity

Hydroxide-catalyzed exchange rate constants were determined for those amides of FK506-binding protein (FKBP12), ubiquitin, and chymotrypsin inhibitor 2 (CI2) that are solvent-accessible in the high-resolution X-ray structures. When combined with previous hydrogen exchange results for the rubredoxin from Pyrococcus furiosus, the acidity of these amides was calculated by continuum dielectric methods as a function of the nonpolarizable electrostatic parameter set, internal dielectric, and the charge distribution of the peptide anion. The CHARMM22 parameter set with an internal dielectric value of 3 and an ab initio-derived anion charge distribution yielded an rmsd value of 7 for the 56 amide exchange rate constants ranging from 10(0.67) to 10(9.0) M(-1) s(-1). The OPLS-AA parameter set yielded comparably robust predictions, while that of PARSE, AMBER parm99, and AMBER ff03 performed more poorly. The small value for the optimal internal dielectric, combined with the brief lifetime of the peptide anion intermediate and the uniformity of the correlation between predicted and observed amide acidities, is consistent with electronic polarizability providing the dominant contribution to dielectric shielding. By construction, nonpolarizable force fields do not model electric field attenuation by electronic polarizability. Accurate prediction of the total electrostatic energy by such force fields necessitates the hyperpolarization of the atomic charge values in order to match the average electric field energy density (1/2)epsilon(tau)E(2)(tau) when epsilon(tau) is set to the in vacuo dielectric value of 1. The resulting predictions of the experimental hydrogen exchange data demonstrate the substantial systematic errors in the predicted electrostatic potential that can arise when dielectric shielding due to electronic polarizability is neglected.