QM/MM Simulations with the Gaussian Electrostatic Model: A Density-based Polarizable Potential.

The use of advanced polarizable potentials in quantum mechanical/molecular mechanical (QM/MM) simulations has been shown to improve the overall accuracy of the calculation. We have developed a density-based potential called the Gaussian electrostatic model (GEM), which has been shown to provide very accurate environments for QM wave functions in QM/MM. In this contribution we present a new implementation of QM/GEM that extends our implementation to include all components (Coulomb, exchange-repulsion, polarization, and dispersion) for the total intermolecular interaction energy in QM/MM calculations, except for the charge-transfer term. The accuracy of the method is tested using a subset of water dimers from the water dimer potential energy surface reported by Babin et al. ( J. Chem. Theory Comput. 2013 9, 5395-5403). Additionally, results of the new implementation are contrasted with results obtained with the classical AMOEBA potential. Our results indicate that GEM provides an accurate MM environment with average root-mean-square error <0.15 kcal/mol for every intermolecular interaction energy component compared with SAPT2+3/aug-cc-pVTZ reference calculations.

S/G-1: an ab initio force-field blending frozen Hermite Gaussian densities and distributed multipoles. Proof of concept and first applications to metal cations.

We demonstrate as a proof of principle the capabilities of a novel hybrid MM'/MM polarizable force field to integrate short-range quantum effects in molecular mechanics (MM) through the use of Gaussian electrostatics. This lead to a further gain in accuracy in the representation of the first coordination shell of metal ions. It uses advanced electrostatics and couples two point dipole polarizable force fields, namely, the Gaussian electrostatic model (GEM), a model based on density fitting, which uses fitted electronic densities to evaluate nonbonded interactions, and SIBFA (sum of interactions between fragments ab initio computed), which resorts to distributed multipoles. To understand the benefits of the use of Gaussian electrostatics, we evaluate first the accuracy of GEM, which is a pure density-based Gaussian electrostatics model on a test Ca(II)-H2O complex. GEM is shown to further improve the agreement of MM polarization with ab initio reference results. Indeed, GEM introduces nonclassical effects by modeling the short-range quantum behavior of electric fields and therefore enables a straightforward (and selective) inclusion of the sole overlap-dependent exchange-polarization repulsive contribution by means of a Gaussian damping function acting on the GEM fields. The S/G-1 scheme is then introduced. Upon limiting the use of Gaussian electrostatics to metal centers only, it is shown to be able to capture the dominant quantum effects at play on the metal coordination sphere. S/G-1 is able to accurately reproduce ab initio total interaction energies within closed-shell metal complexes regarding each individual contribution including the separate contributions of induction, polarization, and charge-transfer. Applications of the method are provided for various systems including the HIV-1 NCp7-Zn(II) metalloprotein. S/G-1 is then extended to heavy metal complexes. Tested on Hg(II) water complexes, S/G-1 is shown to accurately model polarization up to quadrupolar response level. This opens up the possibility of embodying explicit scalar relativistic effects in molecular mechanics thanks to the direct transferability of ab initio pseudopotentials. Therefore, incorporating GEM-like electron density for a metal cation enable the introduction of nonambiguous short-range quantum effects within any point-dipole based polarizable force field without the need of an extensive parametrization.

A finite field method for calculating molecular polarizability tensors for arbitrary multipole rank.

A finite field method for calculating spherical tensor molecular polarizability tensors α(lm;l'm') = ∂Δ(lm)/∂ϕ(l'm')* by numerical derivatives of induced molecular multipole Δ(lm) with respect to gradients of electrostatic potential ϕ(l'm')* is described for arbitrary multipole ranks l and l'. Interconversion formulae for transforming multipole moments and polarizability tensors between spherical and traceless Cartesian tensor conventions are derived. As an example, molecular polarizability tensors up to the hexadecapole-hexadecapole level are calculated for water using the following ab initio methods: Hartree-Fock (HF), Becke three-parameter Lee-Yang-Parr exchange-correlation functional (B3LYP), Møller-Plesset perturbation theory up to second order (MP2), and Coupled Cluster theory with single and double excitations (CCSD). In addition, intermolecular electrostatic and polarization energies calculated by molecular multipoles and polarizability tensors are compared with ab initio reference values calculated by the Reduced Variation Space method for several randomly oriented small molecule dimers separated by a large distance. It is discussed how higher order molecular polarizability tensors can be used as a tool for testing and developing new polarization models for future force fields.

Atomic forces for geometry-dependent point multipole and gaussian multipole models.

In standard treatments of atomic multipole models, interaction energies, total molecular forces, and total molecular torques are given for multipolar interactions between rigid molecules. However, if the molecules are assumed to be flexible, two additional multipolar atomic forces arise because of (1) the transfer of torque between neighboring atoms and (2) the dependence of multipole moment on internal geometry (bond lengths, bond angles, etc.) for geometry-dependent multipole models. In this study, atomic force expressions for geometry-dependent multipoles are presented for use in simulations of flexible molecules. The atomic forces are derived by first proposing a new general expression for Wigner function derivatives partial derivative D(m'm)(l)/partial derivative Omega. The force equations can be applied to electrostatic models based on atomic point multipoles or gaussian multipole charge density. Hydrogen-bonded dimers are used to test the intermolecular electrostatic energies and atomic forces calculated by geometry-dependent multipoles fit to the ab initio electrostatic potential. The electrostatic energies and forces are compared with their reference ab initio values. It is shown that both static and geometry-dependent multipole models are able to reproduce total molecular forces and torques with respect to ab initio, whereas geometry-dependent multipoles are needed to reproduce ab initio atomic forces. The expressions for atomic force can be used in simulations of flexible molecules with atomic multipoles. In addition, the results presented in this work should lead to further development of next generation force fields composed of geometry-dependent multipole models.

Gaussian Multipole Model (GMM).

An electrostatic model based on charge density is proposed as a model for future force fields. The model is composed of a nucleus and a single Slater-type contracted Gaussian multipole charge density on each atom. The Gaussian multipoles are fit to the electrostatic potential (ESP) calculated at the B3LYP/6-31G* and HF/aug-cc-pVTZ levels of theory and tested by comparing electrostatic dimer energies, inter-molecular density overlap integrals, and permanent molecular multipole moments with their respective ab initio values. For the case of water, the atomic Gaussian multipole moments Q(lm) are shown to be a smooth function of internal geometry (bond length and bond angle), which can be approximated by a truncated linear Taylor series. In addition, results are given when the Gaussian multipole charge density is applied to a model for exchange-repulsion energy based on the inter-molecular density overlap.

Conformations and free energy landscapes of polyproline peptides.

The structure of the proline amino acid allows folded polyproline peptides to exist as both left- (PPII) and right-handed (PPI) helices. We have characterized the free energy landscapes of hexamer, nanomer, and tridecamer polyproline peptides in gas phase and implicit water as well as explicit hexane and 1-propanol for the nanomer. To enhance the sampling provided by regular molecular dynamics, we used the recently developed adaptively biased molecular dynamics method, which describes Landau free energy maps in terms of relevant collective variables. These maps, as a function of the collective variables of handedness, radius of gyration, and three others based on the peptide torsion angle omega, were used to determine the relative stability of the different structures, along with an estimate of the transition pathways connecting the different minima. Results show the existence of several metastable isomers and therefore provide a complementary view to experimental conclusions based on photo-induced electron transfer experiments with regard to the existence of stable heterogeneous subpopulations in PPII polyproline.

Structural characterization of the conformational change in calbindin-D28k upon calcium binding using differential surface modification analyzed by mass spectrometry.

Calbindin-D28k is a calcium binding protein with six EF hand domains. Calbindin-D28k is unique in that it functions as both a calcium buffer and a sensor protein. It is found in many tissues, including brain, pancreas, kidney, and intestine, playing important roles in each. Calbindin-D28k is known to bind four calcium ions and upon calcium binding undergoes a conformational change. The structure of apo calbindin-D28k is in an ordered state, transitioning into a disordered state as calcium is bound. Once fully loaded with four calcium ions, it again takes on an ordered state. The solution structure of disulfide-reduced holo-calbindin-D28k has been determined by NMR, while the structure of apo calbindin-D28k has yet to be determined. Differential surface modification of lysine and histidine residues analyzed by mass spectrometry has been used in this study to identify, for the first time, the specific regions of calbindin-D28k undergoing conformational changes between the holo and apo states. Using differential surface modification in combination with mass spectrometry, EF hands 1 and 4 as well as the linkers before EF hand 1 and the linkers between EF hands 4 and 5 and EF hands 5 and 6 were identified as regions of conformational change between apo and holo calbindin-D28k. Under the experimental conditions employed, EF hands 2 and 6, which are known not to bind calcium, were unaffected in either form. EF hand 2 is highly accessible; however, EF hand 6 was determined not to be surface accessible in either form. Previous research has identified a disulfide bond between cysteines 94 and 100 in the holo state. Until now, it was unknown whether this bond also exists in the apo form. Our data confirm the presence of the disulfide bond between cysteines 94 and 100 in the holo form and indicate that there is predominantly no disulfide bond between these residues in the apoprotein.

Reaction mechanism of the epsilon subunit of E. coli DNA polymerase III: insights into active site metal coordination and catalytically significant residues.

The 28 kDa epsilon subunit of Escherichia coli DNA polymerase III is the exonucleotidic proofreader responsible for editing polymerase insertion errors. Here, we study the mechanism by which epsilon carries out the exonuclease activity. We performed quantum mechanics/molecular mechanics calculations on the N-terminal domain containing the exonuclease activity. Both the free-epsilon and a complex epsilon bound to a theta homologue (HOT) were studied. For the epsilon-HOT complex Mg(2+) or Mn(2+) were investigated as the essential divalent metal cofactors, while only Mg(2+) was used for free-epsilon. In all calculations a water molecule bound to the catalytic metal acts as the nucleophile for hydrolysis of the phosphate bond. Initially, a direct proton transfer to H162 is observed. Subsequently, the nucleophilic attack takes place followed by a second proton transfer to E14. Our results show that the reaction catalyzed with Mn(2+) is faster than that with Mg(2+), in agreement with experiment. In addition, the epsilon-HOT complex shows a slightly lower energy barrier compared to free-epsilon. In all cases the catalytic metal is observed to be pentacoordinated. Charge and frontier orbital analyses suggest that charge transfer may stabilize the pentacoordination. Energy decomposition analysis to study the contribution of each residue to catalysis suggests that there are several important residues. Among these, H98, D103, D129, and D146 have been implicated in catalysis by mutagenesis studies. Some of these residues were found to be structurally conserved on human TREX1, the exonuclease domains from E. coli DNA-Pol I, and the DNA polymerase of bacteriophage RB69.

Staggered Mesh Ewald: An extension of the Smooth Particle-Mesh Ewald method adding great versatility.

We draw on an old technique for improving the accuracy of mesh-based field calculations to extend the popular Smooth Particle Mesh Ewald (SPME) algorithm as the Staggered Mesh Ewald (StME) algorithm. StME improves the accuracy of computed forces by up to 1.2 orders of magnitude and also reduces the drift in system momentum inherent in the SPME method by averaging the results of two separate reciprocal space calculations. StME can use charge mesh spacings roughly 1.5× larger than SPME to obtain comparable levels of accuracy; the one mesh in an SPME calculation can therefore be replaced with two separate meshes, each less than one third of the original size. Coarsening the charge mesh can be balanced with reductions in the direct space cutoff to optimize performance: the efficiency of StME rivals or exceeds that of SPME calculations with similarly optimized parameters. StME may also offer advantages for parallel molecular dynamics simulations because it permits the use of coarser meshes without requiring higher orders of charge interpolation and also because the two reciprocal space calculations can be run independently if that is most suitable for the machine architecture. We are planning other improvements to the standard SPME algorithm, and anticipate that StME will work synergistically will all of them to dramatically improve the efficiency and parallel scaling of molecular simulations.

NMDA receptor activation by HIV-Tat protein is clade dependent.

In countries infected with HIV clade B, some patients develop a rapidly progressive dementia that if untreated results in death. In regions of the world infected with HIV clade C, only milder forms of cognitive impairment have been recognized. HIV-infected macrophages are the principal mediators of dementia. HIV clade C, however, efficiently infects macrophages and HIV-infected macrophages are found in the brains of clade C-infected patients. HIV-infected macrophages release Tat protein, which may act directly on neurons to cause toxicity. We found that Tat released from Tat-expressing cells was at least 1000-fold more toxic than recombinant Tat protein. We determined whether Tat could interact with NMDA receptors and whether these interactions are clade dependent. It is demonstrated that Tat binds directly to the NMDA receptor leading to excitotoxicity. The Cys 30-Cys 31 motif in Tat is critical for exciting the NMDA receptor and the Cys31Ser mutation found in clade C Tat has a significantly attenuated neurotoxic response. Through molecular modeling and site-directed mutagenesis, we predict that Cys 31 disrupts the disulfide bond between Cys 744 and Cys 798 on the NR1 subunit of the NMDA receptor by directly interacting with Cys 744 leading to a free thiol group on Cys 798 and subsequent persistent activation of the NMDA receptor.

Identification of a specific motif of the DSS1 protein required for proteasome interaction and p53 protein degradation.

Deleted in Split hand/Split foot 1 (DSS1) was previously identified as a novel 12-O-tetradecanoylphorbol-13-acetate (TPA)-inducible gene with possible involvement in early event of mouse skin carcinogenesis. The mechanisms by which human DSS1 (HsDSS1) exerts its biological effects via regulation of the ubiquitin-proteasome system (UPS) are currently unknown. Here, we demonstrated that HsDSS1 regulates the human proteasome by associating with it in the cytosol and nucleus via the RPN3/S3 subunit of the 19S regulatory particle (RP). Molecular anatomy of HsDSS1 revealed an RPN3/S3-interacting motif (R3IM), located at amino acid residues 15 to 21 of the NH(2) terminus. Importantly, negative charges of the R3IM motif were demonstrated to be required for proteasome interaction and binding to poly-ubiquitinated substrates. Indeed, the R3IM motif of HsDSS1 protein alone was sufficient to replace the ability of intact HsDSS1 protein to pull down proteasome complexes and protein substrates with high-molecular mass ubiquitin conjugates. Interestingly, this interaction is highly conserved throughout evolution from humans to nematodes. Functional study, lowering the levels of the endogenous HsDSS1 using siRNA, indicates that the R3IM/proteasome complex binds and targets p53 for ubiquitin-mediated degradation via gankyrin-MDM2/HDM2 pathway. Most significantly, this work indicates that the R3IM motif of HsDSS1, in conjunction with the complexes of 19S RP and 20S core particle (CP), regulates proteasome interaction through RPN3/S3 molecule, and utilizes a specific subset of poly-ubiquitinated p53 as a substrate.

Calculation of protein-ligand binding free energy by using a polarizable potential.

The binding of charged ligands benzamidine and diazamidine to trypsin was investigated by using a polarizable potential energy function and explicit-water molecular dynamics simulations. The binding free energies were computed from the difference between the free energies of decoupling the ligand from water and protein environments. Both the absolute and the relative free energies from the perturbation simulations agree with experimental measurements to within 0.5 kcal.mol(-1). Comparison of free-energy components sampled from different thermodynamic paths indicates that electrostatics is the main driving force behind benzamidine recognition of trypsin. The contribution of electronic polarization to binding appears to be crucial. By computing the free-energy contribution caused by the polarization between the ligand and its surroundings, we found that polarization has the opposite effect in dissimilar environments. Although polarization favors ligand solvation in water, it weakens the protein-ligand attraction by screening the electrostatic interaction between trypsin and benzamidine. We also examined the relative binding free energies of a benzamidine analog diazamidine to trypsin. The changes in free energy on benzamidine-diazamidine substitution were tens of kilocalories in both water and trypsin environments; however, the change in the total binding free energy is <2 kcal.mol(-1) because of cancellation, consistent with the experimental results. Overall, our results suggest that the use of a polarizable force field, given adequate sampling, is capable of achieving chemical accuracy in molecular simulations of protein-ligand recognition.

Numerical fitting of molecular properties to Hermite Gaussians.

A procedure is presented to fit gridded molecular properties to auxiliary basis sets (ABSs) of Hermite Gaussians, analogous to the density fitting (DF) method (Dunlap; et al. J. Chem. Phys. 1979, 71, 4993). In this procedure, the ab initio calculated properties (density, electrostatic potential, and/or electric field) are fitted via a linear- or nonlinear-least-squares procedure to auxiliary basis sets (ABS). The calculated fitting coefficients from the numerical grids are shown to be more robust than analytic density fitting due to the neglect of the core contributions. The fitting coefficients are tested by calculating intermolecular Coulomb and exchange interactions for a set of dimers. It is shown that the numerical instabilities observed in DF are caused by the attempt of the ABS to fit the core contributions. In addition, this new approach allows us to reduce the number of functions required to obtain an accurate fit. This results in decreased computational cost, which is shown by calculating the Coulomb energy of a 4096 water box in periodic boundary conditions. Using atom centered Hermite Gaussians, this calculation is only 1 order of magnitude slower than conventional atom-centered point charges.

The electrostatic surface term: (I) periodic systems.

The authors propose a new approach to understand the electrostatic surface contributions to the interactions of large but finite periodic distributions of charges. They present a simple method to derive and interpret the surface contribution to any electrostatic field produced by a periodic distribution of charges. They discuss the physical and mathematical interpretations of this term. They present several examples and physical details associated with the calculation of the surface term. Finally, they provide a simple derivation of the surface contribution to the virial. This term does not disappear even if tinfoil boundary conditions are applied.

Anisotropic, Polarizable Molecular Mechanics Studies of Inter- and Intramolecular Interactions and Ligand-Macromolecule Complexes. A Bottom-Up Strategy.

We present an overview of the SIBFA polarizable molecular mechanics procedure, which is formulated and calibrated on the basis of quantum chemistry (QC). It embodies nonclassical effects such as electrostatic penetration, exchange-polarization, and charge transfer. We address the issues of anisotropy, nonadditivity, and transferability by performing parallel QC computations on multimolecular complexes. These encompass multiply H-bonded complexes and polycoordinated complexes of divalent cations. Recent applications to the docking of inhibitors to Zn-metalloproteins are presented next, namely metallo-beta-lactamase, phosphomannoisomerase, and the nucleocapsid of the HIV-1 retrovirus. Finally, toward third-generation intermolecular potentials based on density fitting, we present the development of a novel methodology, the Gaussian electrostatic model (GEM), which relies on ab initio-derived fragment electron densities to compute the components of the total interaction energy. As GEM offers the possibility of a continuous electrostatic model going from distributed multipoles to densities, it allows an inclusion of short-range quantum effects in the molecular mechanics energies. The perspectives of an integrated SIBFA/GEM/QM procedure are discussed.

Mapping the dimer interface in the C-terminal domains of the yeast MLH1-PMS1 heterodimer.

Yeast MutLalpha is a heterodimer of MLH1 and PMS1 that participates in a variety of DNA transactions, including DNA mismatch repair. Formation of the MutLalpha heterodimer requires that the C-terminal domains of MLH1 and PMS1 interact in a manner that is not yet fully understood. Here we investigate the interactions involved in heterodimerization. Using protein surface modification and mass spectrometry, we identify numerous lysine residues that are exposed to solvent in monomeric MLH1. A corresponding analysis of the MLH1-PMS1 heterodimer reveals that three of these exposed residues, K665, K675, and K704, are no longer solvent accessible in the heterodimer, suggesting that they are within the dimer interface. We refine secondary structure predictions and sequence alignments of C-terminal residues of seven eukaryotic MutL homologues and then develop homology models for the N- and C-terminal domains of MLH1. On the basis of this information, we present a model for interaction of the C-terminal domains of MLH1 and PMS1.

The free energy landscape of small peptides as obtained from metadynamics with umbrella sampling corrections.

There is considerable interest in developing methodologies for the accurate evaluation of free energies, especially in the context of biomolecular simulations. Here, we report on a reexamination of the recently developed metadynamics method, which is explicitly designed to probe "rare events" and areas of phase space that are typically difficult to access with a molecular dynamics simulation. Specifically, we show that the accuracy of the free energy landscape calculated with the metadynamics method may be considerably improved when combined with umbrella sampling techniques. As test cases, we have studied the folding free energy landscape of two prototypical peptides: Ace-(Gly)(2)-Pro-(Gly)(3)-Nme in vacuo and trialanine solvated by both implicit and explicit water. The method has been implemented in the classical biomolecular code AMBER and is to be distributed in the next scheduled release of the code.

Generalization of the Gaussian electrostatic model: extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods.

The simulation of biological systems by means of current empirical force fields presents shortcomings due to their lack of accuracy, especially in the description of the nonbonded terms. We have previously introduced a force field based on density fitting termed the Gaussian electrostatic model-0 (GEM-0) J.-P. Piquemal et al. [J. Chem. Phys. 124, 104101 (2006)] that improves the description of the nonbonded interactions. GEM-0 relies on density fitting methodology to reproduce each contribution of the constrained space orbital variation (CSOV) energy decomposition scheme, by expanding the electronic density of the molecule in s-type Gaussian functions centered at specific sites. In the present contribution we extend the Coulomb and exchange components of the force field to auxiliary basis sets of arbitrary angular momentum. Since the basis functions with higher angular momentum have directionality, a reference molecular frame (local frame) formalism is employed for the rotation of the fitted expansion coefficients. In all cases the intermolecular interaction energies are calculated by means of Hermite Gaussian functions using the McMurchie-Davidson [J. Comput. Phys. 26, 218 (1978)] recursion to calculate all the required integrals. Furthermore, the use of Hermite Gaussian functions allows a point multipole decomposition determination at each expansion site. Additionally, the issue of computational speed is investigated by reciprocal space based formalisms which include the particle mesh Ewald (PME) and fast Fourier-Poisson (FFP) methods. Frozen-core (Coulomb and exchange-repulsion) intermolecular interaction results for ten stationary points on the water dimer potential-energy surface, as well as a one-dimensional surface scan for the canonical water dimer, formamide, stacked benzene, and benzene water dimers, are presented. All results show reasonable agreement with the corresponding CSOV calculated reference contributions, around 0.1 and 0.15 kcal/mol error for Coulomb and exchange, respectively. Timing results for single Coulomb energy-force calculations for (H(2)O)(n), n=64, 128, 256, 512, and 1024, in periodic boundary conditions with PME and FFP at two different rms force tolerances are also presented. For the small and intermediate auxiliaries, PME shows faster times than FFP at both accuracies and the advantage of PME widens at higher accuracy, while for the largest auxiliary, the opposite occurs.

Towards accurate solvation dynamics of divalent cations in water using the polarizable amoeba force field: From energetics to structure.

Molecular dynamics simulations were performed using a modified amoeba force field to determine hydration and dynamical properties of the divalent cations Ca2+ and Mg2+. The extension of amoeba to divalent cations required the introduction of a cation specific parametrization. To accomplish this, the Thole polarization damping model parametrization was modified based on the ab initio polarization energy computed by a constrained space orbital variation energy decomposition scheme. Excellent agreement has been found with condensed phase experimental results using parameters derived from gas phase ab initio calculations. Additionally, we have observed that the coordination of the calcium cation is influenced by the size of the periodic water box, a recurrent issue in first principles molecular dynamics studies.

Simulation of Ca2+ and Mg2+ solvation using polarizable atomic multipole potential.

The alkaline earth metals calcium and magnesium are critically involved in many biomolecular processes. To understand the hydration thermodynamics of these ions, we have performed molecular dynamics simulations using a polarizable potential. Particle-mesh Ewald for point multipoles has been applied to the calculation of electrostatic interactions. The parameters in this model have been determined from an ab initio quantum mechanical calculation of dimer interactions between ions and water. Two methods for ion solvation free energy calculation, free energy perturbation, and the Bennett acceptance ratio have been compared. Both predict results consistent with other theoretical estimations while the Bennett approach leads to a much smaller statistical error. Based on the Born theory and the ion-oxygen radial distribution functions, we estimate the effective size of the ions in solution, concluding that K(+) > Na(+) congruent with Ca(2+) > Mg(2+). There appears to be much stronger perturbation in water structure, dynamics, and dipole moment around the divalent cations than the monovalent K(+) and Na(+). The average water coordination numbers for Ca(2+) and Mg(2+) are 7.3 and 6, respectively. The lifetime of water molecules in the first solvation shell of Mg(2+) is on the order of hundreds of picoseconds, in contrast to only few picoseconds for Ca(2+), K(+), or Na(+).