Current computer simulations of biomolecules typically make use of classical molecular dynamics methods, as a very large number (tens to hundreds of thousands) of atoms are involved over timescales of many nanoseconds. The methodology for treating short-range bonded and van der Waals interactions has matured. However, long-range electrostatic interactions still represent a bottleneck in simulations. In this article, we introduce the basic issues for an accurate representation of the relevant electrostatic interactions. In spite of the huge computational time demanded by most biomolecular systems, it is no longer necessary to resort to uncontrolled approximations such as the use of cutoffs. In particular, we discuss the Ewald summation methods, the fast particle mesh methods, and the fast multipole methods. We also review recent efforts to understand the role of boundary conditions in systems with long-range interactions, and conclude with a short perspective on future trends.

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-β-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.

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(+).

LI3@niehs.nih.gov

Total intermolecular interaction energies are determined with a first version of the Gaussian electrostatic model (GEM-0), a force field based on a density fitting approach using s-type Gaussian functions. The total interaction energy is computed in the spirit of the sum of interacting fragment ab initio (SIBFA) force field by separately evaluating each one of its components: electrostatic (Coulomb), exchange repulsion, polarization, and charge transfer intermolecular interaction energies, in order to reproduce reference constrained space orbital variation (CSOV) energy decomposition calculations at the B3LYP/aug-cc-pVTZ level. The use of an auxiliary basis set restricted to spherical Gaussian functions facilitates the rotation of the fitted densities of rigid fragments and enables a fast and accurate density fitting evaluation of Coulomb and exchangerepulsion energy, the latter using the overlap model introduced by Wheatley and Price [Mol. Phys. 69, 50718 (1990)]. The SIBFA energy scheme for polarization and charge transfer has been implemented using the electric fields and electrostatic potentials generated by the fitted densities. GEM-0 has been tested on ten stationary points of the water dimer potential energy surface and on three water clusters (n=16,20,64). The results show very good agreement with density functional theory calculations, reproducing the individual CSOV energy contributions for a given interaction as well as the B3LYP total interaction energies with errors below k B T at room temperature. Preliminary results for Coulomb and exchange-repulsion energies of metal cation complexes and coupled cluster singles doubles electron densities are discussed.

The accurate simulation of biologically active macromolecules faces serious limitations that originate in the treatment of electrostatics in the empirical force fields. The current use of "partial charges" is a significant source of errors, since these vary widely with different conformations. By contrast, the molecular electrostatic potential (MEP) obtained through the use of a distributed multipole moment description, has been shown to converge to the quantum MEP outside the van der Waals surface, when higher order multipoles are used. However, in spite of the considerable improvement to the representation of the electronic cloud, higher order multipoles are not part of current classical biomolecular force fields due to the excessive computational cost. In this paper we present an efficient formalism for the treatment of higher order multipoles in Cartesian tensor formalism. The Ewald "direct sum" is evaluated through a McMurchie-Davidson formalism [L. McMurchie and E. Davidson, J. Comput. Phys. 26, 218 (1978)]. The "reciprocal sum" has been implemented in three different ways: using an Ewald scheme, a particle mesh Ewald (PME) method, and a multigrid-based approach. We find that even though the use of the McMurchie-Davidson formalism considerably reduces the cost of the calculation with respect to the standard matrix implementation of multipole interactions, the calculation in direct space remains expensive. When most of the calculation is moved to reciprocal space via the PME method, the cost of a calculation where all multipolar interactions (up to hexadecapole-hexadecapole) are included is only about 8.5 times more expensive than a regular AMBER 7 [D. A. Pearlman et al., Comput. Phys. Commun. 91, 1 (1995)] implementation with only charge-charge interactions. The multigrid implementation is slower but shows very promising results for parallelization. It provides a natural way to interface with continuous, Gaussian-based electrostatics in the future. It is hoped that this new formalism will facilitate the systematic implementation of higher order multipoles in classical biomolecular force fields.

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.simulation ͉ molecular dynamics ͉ trypsin ͉ benzamidine ͉ force field S pecific recognition of ligands by proteins is at the core of many crucial biological functions and systems such as enzyme catalysis and intracellular signaling. Binding affinity characterizes the strength of such recognition. With the recent advancements in computing, prediction of the binding affinity based on physical principles of molecular interaction has come to the forefront of active research and has been the subject of regular reviews (1-5). All-atom molecular dynamics (MD) simulation with explicit solvent, coupled with efficient free-energy sampling algorithms, can potentially offer accurate prediction of binding free energies of ligands to proteins (5). Common free-energy simulation algorithms include the double-decoupling method (DDM) and potential of mean force approach (PMF). Free-energy perturbation (FEP), thermodynamic integration (TI), or umbrella sampling can be used to compute free-energy differences in either DDM or PMF. It has been argued that DDM is problematic for charged systems, because the binding free energy is computed as a small difference between two large solvation energies in water and in protein (6). However, the PMF approach does not quantify absolute solvation energies of ligand, which makes it difficult to detect potential problems in treatment of long-range effect and bou...

Little is known about the mechanisms that regulate differential transactivation by p53. We developed a system in the yeast Saccharomyces cerevisiae that addresses p53 transactivation capacity from 26 different p53 response elements (REs) under conditions where all other factors, such as chromatin, are kept constant. The system relies on a tightly regulated promoter (rheostatable) that can provide for a broad range of p53 expression. The p53 transactivation capacity toward each 20-to 22-bp-long RE could be ranked by using a simple phenotypic assay. Surprisingly, there was as much as a 1,000-fold difference in transactivation. There was no correlation between the functional rank and statistical predictions of binding energy of the REs. Instead we found that the central sequence element in an RE greatly affects p53 transactivation capacity, possibly because of DNA structural properties. Our results suggest that intrinsic DNA binding affinity and p53 protein levels are important contributors to p53-induced differential transactivation. These results are also relevant to understanding the regulation by other families of transcription factors that recognize several sequence-related response elements and/or have tightly regulated expression. We found that p53 had weak activity towards half the apoptotic REs. In addition, p53 alleles associated with familial breast cancer, previously classified as wild type, showed subtle differences in transactivation capacity towards several REs.