Ewald-based methods for Gaussian integral evaluation: application to a new parameterization of GEM.

The development of accurate potentials for computational simulations has been an active area of research. Our group has been involved in the development of the Gaussian electrostatic model (GEM), a force field based on molecular densities. The philosophy of GEM is based on the pioneering work of N. Gresh and co-workers of the reproduction of individual inter-molecular interaction components obtained from quantum mechanical (QM) energy decomposition analysis (EDA). The molecular densities used in GEM are represented by fitting accurate QM molecular densities using auxiliary basis sets (comprised of Hermite Gaussians). The use of these molecular densities results in the need to evaluate a large number of Gaussian integrals. We have previously shown that the particle-mesh Ewald (PME), and fast Fourier Poisson (FFP) methods can be used for efficiently evaluating these types of integrals. Here, we present the latest parameterization of GEM* and its application for an extensive study of PME and FFP for molecular dynamics (MD) simulations using a hybrid version of our potential, GEM*. The temperature dependence of various bulk properties is presented and discussed, as well as the effect of various parameters affecting the performance/accuracy of both methods.

Long-range electrostatic corrections in multipolar/polarizable QM/MM simulations.

Taking long-range electrostatic effects into account in classical and hybrid quantum mechanics-molecular mechanics (QM/MM) simulations is necessary for an accurate description of the system under study. We have recently developed a method, termed long-range electrostatic corrections (LREC), for monopolar QM/MM calculations. Here, we present an extension of LREC for multipolar/polarizable QM/MM simulations within the LICHEM software package. Reaction barriers and QM-MM interaction energies converge with a LREC cutoff between 20 and 25 Å, in agreement with our previous results. Additionally, the LREC approach for the QM-MM interactions can be smoothly combined with standard shifting or Ewald summation methods in the MM calculations. We recommend the use of QM(LREC)/MM(PME), where the QM region is treated with LREC and the MM region is treated with particle mesh Ewald (PME) summation. This combination is an excellent compromise between simplicity, speed, and accuracy for large QM/MM simulations.

A new smoothing function to introduce long-range electrostatic effects in QM/MM calculations.

A new method to account for long range electrostatic contributions is proposed and implemented for quantum mechanics/molecular mechanics long range electrostatic correction (QM/MM-LREC) calculations. This method involves the use of the minimum image convention under periodic boundary conditions and a new smoothing function for energies and forces at the cutoff boundary for the Coulomb interactions. Compared to conventional QM/MM calculations without long-range electrostatic corrections, the new method effectively includes effects on the MM environment in the primary image from its replicas in the neighborhood. QM/MM-LREC offers three useful features including the avoidance of calculations in reciprocal space (k-space), with the concomitant avoidance of having to reproduce (analytically or approximately) the QM charge density in k-space, and the straightforward availability of analytical Hessians. The new method is tested and compared with results from smooth particle mesh Ewald (PME) for three systems including a box of neat water, a double proton transfer reaction, and the geometry optimization of the critical point structures for the rate limiting step of the DNA dealkylase AlkB. As with other smoothing or shifting functions, relatively large cutoffs are necessary to achieve comparable accuracy with PME. For the double-proton transfer reaction, the use of a 22 Å cutoff shows a close reaction energy profile and geometries of stationary structures with QM/MM-LREC compared to conventional QM/MM with no truncation. Geometry optimization of stationary structures for the hydrogen abstraction step by AlkB shows some differences between QM/MM-LREC and the conventional QM/MM. These differences underscore the necessity of the inclusion of the long-range electrostatic contribution.

GEM*: A Molecular Electronic Density-Based Force Field for Molecular Dynamics Simulations.

GEM*, a force field that combines Coulomb and Exchange terms calculated with Hermite Gaussians with the polarization, bonded, and modified van der Waals terms from AMOEBA is presented. GEM* is tested on an initial water model fitted at the same level as AMOEBA. The integrals required for the evaluation of the intermolecular Coulomb interactions are efficiently evaluated by means of reciprocal space methods. The GEM* water model is tested by comparing energies and forces for a series of water oligomers and MD simulations. Timings for GEM* compared to AMOEBA are presented and discussed.

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.

Simulations of a protein crystal with a high resolution X-ray structure: evaluation of force fields and water models.

We use classical molecular dynamics and 16 combinations of force fields and water models to simulate a protein crystal observed by room-temperature X-ray diffraction. The high resolution of the diffraction data (0.96 Å) and the simplicity of the crystallization solution (nearly pure water) make it possible to attribute any inconsistencies between the crystal structure and our simulations to artifacts of the models rather than inadequate representation of the crystal environment or uncertainty in the experiment. All simulations were extended for 100 ns of production dynamics, permitting some long-time scale artifacts of each model to emerge. The most noticeable effect of these artifacts is a model-dependent drift in the unit cell dimensions, which can become as large as 5% in certain force fields; the underlying cause is the replacement of native crystallographic contacts with non-native ones, which can occur with heterogeneity (loss of crystallographic symmetry) in simulations with some force fields. We find that the AMBER FF99SB force field maintains a lattice structure nearest that seen in the X-ray data, and produces the most realistic atomic fluctuations (by comparison to crystallographic B-factors) of all the models tested. We find that the choice of water model has a minor effect in comparison to the choice of protein model. We also identify a number of artifacts that occur throughout all of the simulations: excessive formation of hydrogen bonds or salt bridges between polar groups and loss of hydrophobic interactions. This study is intended as a foundation for future work that will identify individual parameters in each molecular model that can be modified to improve their representations of protein structure and thermodynamics.

Recent estimates of the structure of the factor VIIa (FVIIa)/tissue factor (TF) and factor Xa (FXa) ternary complex.

The putative structure of the Tissue Factor/Factor VIIa/Factor Xa (TF/FVIIa/FXa) ternary complex is reconsidered. Two independently derived docking models proposed in 2003 (one for our laboratory: CHeA and one from the Scripps laboratory: Ss) are dynamically equilibrated for over 10 ns in an electrically neutral solution using all-atom molecular dynamics. Although the dynamical models (CHeB and Se) differ in atomic detail, there are similarities in that TF is found to interact with the gamma-carboxyglutamic acid (Gla) and Epidermal Growth Factor-like 1 (EGF-1) domains of FXa, and FVIIa is found to interact with the Gla, EGF-2 and serine protease (SP) domains of FXa in both models. FVIIa does not interact with the FXa EGF-1 domain in Se and the EGF domains of FVIIa do not interact with FXa in the CHeB. Both models are consistent with experimentally suggested contacts between the SP domain of FVIIa with the EGF-2 and SP domains of FXa.

Effects of early and late rest breaks during training on overnight memory consolidation of a keyboard melody.

In two experiments, we tested the extent to which overnight procedural memory consolidation is affected by extended rest breaks during training. In the first experiment, nonmusicians practiced a 5-element keypress sequence with their nondominant hand in 12 30-s practice intervals separated by 30-s pauses. In the second experiment, nonpianist musicians practiced a 13-note keyboard melody using the same procedures. In both experiments, approximately one-third of the subjects took a 5-min break after the first three blocks of practice; another third took a break after nine blocks of practice; the remaining participants did not take an extended break. All were trained in the evening and were retested the following morning. Participants in both experiments made dramatic improvements over the course of the training and retest sessions, and participants who took an extended rest break early in practice made the largest gains in performance between the end of training and the beginning of retest.

Trypsin-ligand binding free energies from explicit and implicit solvent simulations with polarizable potential.

We have calculated the binding free energies of a series of benzamidine-like inhibitors to trypsin with a polarizable force field using both explicit and implicit solvent approaches. Free energy perturbation has been performed for the ligands in bulk water and in protein complex with molecular dynamics simulations. The binding free energies calculated from explicit solvent simulations are well within the accuracy of experimental measurement and the direction of change is predicted correctly in all cases. We analyzed the molecular dipole moments of the ligands in gas, water and protein environments. Neither binding affinity nor ligand solvation free energy in bulk water shows much dependence on the molecular dipole moments of the ligands. Substitution of the aromatic or the charged group in the ligand results in considerable change in the solvation energy in bulk water and protein whereas the binding affinity varies insignificantly due to cancellation. The effect of chemical modification on ligand charge distribution is mostly local. Replacing benzene with diazine has minimal impact on the atomic multipoles at the amidinium group. We have also utilized an implicit solvent based end-state approach to evaluate the binding free energies of these inhibitors. In this approach, the polarizable multipole model combined with Poisson-Boltzmann/surface area (PMPB/SA) provides the electrostatic interaction energy and the polar solvation free energy. Overall the relative binding free energies obtained from the MM-PMPB/SA model are in good agreement with the experimental data.

A computational modeling and molecular dynamics study of the Michaelis complex of human protein Z-dependent protease inhibitor (ZPI) and factor Xa (FXa).

Protein Z-dependent protease inhibitor (ZPI) and antithrombin III (AT3) are members of the serpin superfamily of protease inhibitors that inhibit factor Xa (FXa) and other proteases in the coagulation pathway. While experimental structural information is available for the interaction of AT3 with FXa, at present there is no structural data regarding the interaction of ZPI with FXa, and the precise role of this interaction in the blood coagulation pathway is poorly understood. In an effort to gain a structural understanding of this system, we have built a solvent equilibrated three-dimensional structural model of the Michaelis complex of human ZPI/FXa using homology modeling, protein-protein docking and molecular dynamics simulation methods. Preliminary analysis of interactions at the complex interface from our simulations suggests that the interactions of the reactive center loop (RCL) and the exosite surface of ZPI with FXa are similar to those observed from X-ray crystal structure-based simulations of AT3/FXa. However, detailed comparison of our modeled structure of ZPI/FXa with that of AT3/FXa points to differences in interaction specificity at the reactive center and in the stability of the inhibitory complex, due to the presence of a tyrosine residue at the P1 position in ZPI, instead of the P1 arginine residue in AT3. The modeled structure also shows specific structural differences between AT3 and ZPI in the heparin-binding and flexible N-terminal tail regions. Our structural model of ZPI/FXa is also compatible with available experimental information regarding the importance for the inhibitory action of certain basic residues in FXa.

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.

Computational study of the putative active form of protein Z (PZa): sequence design and structural modeling.

Although protein Z (PZ) has a domain arrangement similar to the essential coagulation proteins FVII, FIX, FX, and protein C, its serine protease (SP)-like domain is incomplete and does not exhibit proteolytic activity. We have generated a trial sequence of putative activated protein Z (PZa) by identifying amino acid mutations in the SP-like domain that might reasonably resurrect the serine protease catalytic activity of PZ. The structure of the activated form was then modeled based on the proposed sequence using homology modeling and solvent-equilibrated molecular dynamics simulations. In silico docking of inhibitors of FVIIa and FXa to the putative active site of equilibrated PZa, along with structural comparison with its homologous proteins, suggest that the designed PZa can possibly act as a serine protease.

Vulnerability in Popular Molecular Dynamics Packages Concerning Langevin and Andersen Dynamics.

We report a serious problem associated with a number of current implementations of Andersen and Langevin dynamics algorithms. When long simulations are run in many segments, it is sometimes possible to have a repeating sequence of pseudorandom numbers enter the calcuation. We show that, if the sequence repeats rapidly, the resulting artifacts can quickly denature biomolecules and are then easily detectable. However, if the sequence repeats less frequently, the artifacts become subtle and easily overlooked. We derive a formula for the underlying cause of artifacts in the case of the Langevin thermostat, and find it vanishes slowly as the inverse square root of the number of time steps per simulation segment. Numerous examples of simulation artifacts are presented, including dissociation of a tetrameric protein after 110 ns of dynamics, reductions in atomic fluctuations for a small protein in implicit solvent, altered thermodynamic properties of a box of water molecules, and changes in the transition free energies between dihedral angle conformations. Finally, in the case of strong thermocoupling, we link the observed artifacts to previous work in nonlinear dynamics and show that it is possible to drive a 20-residue, implicitly solvated protein into periodic trajectories if the thermostat is not used properly. Our findings should help other investigators re-evaluate simulations that may have been corrupted and obtain more accurate results.

Cough syncope from constrictive pericarditis: a case report.

Loss of consciousness associated in some people with coughing is likely due to the hemodynamic changes induced by this exaggerated Valsalva manoeuvre. Cough syncope classically occurs in people with obstructive pulmonary disorders and it is unusual to find a specific correctable lesion precipitating it. A previously healthy man who developed tussive syncope in association with constrictive pericarditis and subsequently was cured by pericardiectomy is described. Although cough is a recognized feature of this cardiovascular condition, associated cough syncope has not previously been reported. However, the hemodynamic changes induced by constriction of cardiac filling understandably may further reduce cerebral perfusion during coughing. Not only does this case present a novel cause of this interesting syncopal syndrome, but it supports the basic pathophysiology thought to underlie it.