PUPIL: A Software Integration System for Multi-Scale QM/MM-MD Simulations and Its Application to Biomolecular Systems.

PUPIL (Program for User Package Interfacing and Linking) implements a distinctive multi-scale approach to hybrid quantum mechanical/molecular mechanical molecular dynamics (QM/MM-MD) simulations. Originally developed to interface different external programs for multi-scale simulation with applications in the materials sciences, PUPIL is finding increasing use in the study of complex biological systems. Advanced MD techniques from the external packages can be applied readily to a hybrid QM/MM treatment in which the forces and energy for the QM region can be computed by any of the QM methods available in any of the other external packages. Here, we give a survey of PUPIL design philosophy, main features, and key implementation decisions, with an orientation to biomolecular simulation. We discuss recently implemented features which enable highly realistic simulations of complex biological systems which have more than one active site that must be treated concurrently. Examples are given.

Rational Design of Particle Mesh Ewald Compatible Lennard-Jones Parameters for +2 Metal Cations in Explicit Solvent.

Metal ions play significant roles in biological systems. Accurate molecular dynamics (MD) simulations on these systems require a validated set of parameters. Although there are more detailed ways to model metal ions, the nonbonded model, which employs a 12-6 Lennard-Jones (LJ) term plus an electrostatic potential is still widely used in MD simulations today due to its simple form. However, LJ parameters have limited transferability due to different combining rules, various water models and diverse simulation methods. Recently, simulations employing a Particle Mesh Ewald (PME) treatment for long-range electrostatics have become more and more popular owing to their speed and accuracy. In the present work we have systematically designed LJ parameters for 24 +2 metal (M(II)) cations to reproduce different experimental properties appropriate for the Lorentz-Berthelot combining rules and PME simulations. We began by testing the transferability of currently available M(II) ion LJ parameters. The results showed that there are differences between simulations employing Ewald summation with other simulation methods and that it was necessary to design new parameters specific for PME based simulations. Employing the thermodynamic integration (TI) method and performing periodic boundary MD simulations employing PME, allowed for the systematic investigation of the LJ parameter space. Hydration free energies (HFEs), the ion-oxygen distance in the first solvation shell (IOD) and coordination numbers (CNs) were obtained for various combinations of the parameters of the LJ potential for four widely used water models (TIP3P, SPC/E, TIP4P and TIP4PEW). Results showed that the three simulated properties were highly correlated. Meanwhile, M(II) ions with the same parameters in different water models produce remarkably different HFEs but similar structural properties. It is difficult to reproduce various experimental values simultaneously because the nonbonded model underestimates the interaction between the metal ions and water molecules at short range. Moreover, the extent of underestimation increases successively for the TIP3P, SPC/E, TIP4PEW and TIP4P water models. Nonetheless, we fitted a curve to describe the relationship between ε (the well depth) and radius (Rmin/2) from experimental data on noble gases to facilitate the generation of the best possible compromise models. Hence, by targeting different experimental values, we developed three sets of parameters for M(II) cations for three different water models (TIP3P, SPC/E and TIP4PEW). These parameters we feel represent the best possible compromise that can be achieved using the nonbonded model for the ions in combination with simple water models. From a computational uncertainty analysis we estimate that the uncertainty in our computed HFEs is on the order of ±1kcal/mol. Further improvements will require more advanced non-bonded models likely with inclusion of polarization.

Wide-open flaps are key to urease activity.

Substrate ingress and product egress from the active site of urease is tightly controlled by an active-site flap. Molecular dynamics simulations of urease have revealed a previously unobserved wide-open flap state that, unlike the well-characterized closed and open states, allows ready access to the metal cluster in the active site. This state is easily reached from the open state via low free energy barriers. Additionally, we have found that even when the flap is closed, a region of the binding pocket is solvent-exposed, leading to the hypothesis that it may act as a substrate/product reservoir. The newly identified wide-open state offers further opportunities for small-molecule drug discovery by defining a more extensive active-site pocket than has been previously described.

Comment on "A minimal implementation of the AMBER-GAUSSIAN interface for ab initio QM/MM-MD simulation".

We comment upon the recent critique of use of the Program for User Package Interfacing and Linking (PUPIL) system for linking AMBER and GAUSSIAN in a multiscale quantum mechanical/molecular mechanics (QM/MM) simulation (Okamoto et al., J. Comput. Chem. 2011, 32, 932). Specifically, their method for computing forces on the MM particles from the QM region via the GAUSSIAN-03 electrical field was already implemented in PUPIL version 1.3, publicly available beginning December 2009. Some other doubtful characterizations of PUPIL are discussed briefly in the context of current awareness of open-source codes more generally.

Formal Estimation of Errors in Computed Absolute Interaction Energies of Protein-ligand Complexes.

A largely unsolved problem in computational biochemistry is the accurate prediction of binding affinities of small ligands to protein receptors. We present a detailed analysis of the systematic and random errors present in computational methods through the use of error probability density functions, specifically for computed interaction energies between chemical fragments comprising a protein-ligand complex. An HIV-II protease crystal structure with a bound ligand (indinavir) was chosen as a model protein-ligand complex. The complex was decomposed into twenty-one (21) interacting fragment pairs, which were studied using a number of computational methods. The chemically accurate complete basis set coupled cluster theory (CCSD(T)/CBS) interaction energies were used as reference values to generate our error estimates. In our analysis we observed significant systematic and random errors in most methods, which was surprising especially for parameterized classical and semiempirical quantum mechanical calculations. After propagating these fragment-based error estimates over the entire protein-ligand complex, our total error estimates for many methods are large compared to the experimentally determined free energy of binding. Thus, we conclude that statistical error analysis is a necessary addition to any scoring function attempting to produce reliable binding affinity predictions.

The energy computation paradox and ab initio protein folding.

The routine prediction of three-dimensional protein structure from sequence remains a challenge in computational biochemistry. It has been intuited that calculated energies from physics-based scoring functions are able to distinguish native from nonnative folds based on previous performance with small proteins and that conformational sampling is the fundamental bottleneck to successful folding. We demonstrate that as protein size increases, errors in the computed energies become a significant problem. We show, by using error probability density functions, that physics-based scores contain significant systematic and random errors relative to accurate reference energies. These errors propagate throughout an entire protein and distort its energy landscape to such an extent that modern scoring functions should have little chance of success in finding the free energy minima of large proteins. Nonetheless, by understanding errors in physics-based score functions, they can be reduced in a post-hoc manner, improving accuracy in energy computation and fold discrimination.

Novel acyclic diaminocarbene ligands with increased steric demand and their application in gold catalysis.

Sterically demanding and conformationally stable N,N'-ditertiaryalkyl-N,N'-diphenyl acyclic diaminocarbenes (ADCs) were developed. Bulky ADC-Au catalysts not only showed competitive reactivities in hydroamination and enyne cyclization but also demonstrated unique ligand properties different from bulky N-heterocyclic carbene (NHC) counterparts.