Alexey V. Onufriev

The Amber biomolecular simulation programs.

We describe the development, current features, and some directions for future development of the Amber package of computer programs. This package evolved from a program that was constructed in the late 1970s to do Assisted Model Building with Energy Refinement, and now contains a group of programs embodying a number of powerful tools of modern computational chemistry, focused on molecular dynamics and free energy calculations of proteins, nucleic acids, and carbohydrates.

Exploring protein native states and large‐scale conformational changes with a modified generalized born model

Implicit solvation models provide, for many applications, a reasonably accurate and computationally effective way to describe the electrostatics of aqueous solvation. Here, a popular analytical Generalized Born (GB) solvation model is modified to improve its accuracy in calculating the solvent polarization part of free energy changes in large‐scale conformational transitions, such as protein folding. In contrast to an earlier GB model (implemented in the AMBER‐6 program), the improved version does not overstabilize the native structures relative to the finite‐difference Poisson–Boltzmann continuum treatment. In addition to improving the energy balance between folded and unfolded conformers, the algorithm (available in the AMBER‐7 and NAB molecular modeling packages) is shown to perform well in more than 50 ns of native‐state molecular dynamics (MD) simulations of thioredoxin, protein‐A, and ubiquitin, as well as in a simulation of Barnase/Barstar complex formation. For thioredoxin, various combinations of input parameters have been explored, such as the underlying gas‐phase force fields and the atomic radii. The best performance is achieved with a previously proposed modification to the torsional potential in the Amber ff99 force field, which yields stable native trajectories for all of the tested proteins, with backbone root‐mean‐square deviations from the native structures being ∼1.5 Å after 6 ns of simulation time. The structure of Barnase/Barstar complex is regenerated, starting from an unbound state, to within 1.9 Å relative to the crystal structure of the complex. Proteins 2004;55:000–000.

H++: a server for estimating pKas and adding missing hydrogens to macromolecules.

The structure and function of macromolecules depend critically on the ionization (protonation) states of their acidic and basic groups. A number of existing practical methods predict protonation equilibrium pK constants of macromolecules based upon their atomic resolution Protein Data Bank (PDB) structures; the calculations are often performed within the framework of the continuum electrostatics model. Unfortunately, these methodologies are complex, involve multiple steps and require considerable investment of effort. Our web server provides access to a tool that automates this process, allowing both experts and novices to quickly obtain estimates of pKs as well as other related characteristics of biomolecules such as isoelectric points, titration curves and energies of protonation microstates. Protons are added to the input structure according to the calculated ionization states of its titratable groups at the user-specified pH; the output is in the PQR (PDB + charges + radii) format. In addition, corresponding coordinate and topology files are generated in the format supported by the molecular modeling package AMBER. The server is intended for a broad community of biochemists, molecular modelers, structural biologists and drug designers; it can also be used as an educational tool in biochemistry courses.

Performance comparison of generalized born and Poisson methods in the calculation of electrostatic solvation energies for protein structures

This study compares generalized Born (GB) and Poisson (PB) methods for calculating electrostatic solvation energies of proteins. A large set of GB and PB implementations from our own laboratories as well as others is applied to a series of protein structure test sets for evaluating the performance of these methods. The test sets cover a significant range of native protein structures of varying size, fold topology, and amino acid composition as well as nonnative extended and misfolded structures that may be found during structure prediction and folding/unfolding studies. We find that the methods tested here span a wide range from highly accurate and computationally demanding PB‐based methods to somewhat less accurate but more affordable GB‐based approaches and a few fast, approximate PB solvers. Compared with PB solvation energies, the latest, most accurate GB implementations were found to achieve errors of 1% for relative solvation energies between different proteins and 0.4% between different conformations of the same protein. This compares to accurate PB solvers that produce results with deviations of less than 0.25% between each other for both native and nonnative structures. The performance of the best GB methods is discussed in more detail for the application for force field‐based minimizations or molecular dynamics simulations.

H++ 3.0: automating pK prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulations

The accuracy of atomistic biomolecular modeling and simulation studies depend on the accuracy of the input structures. Preparing these structures for an atomistic modeling task, such as molecular dynamics (MD) simulation, can involve the use of a variety of different tools for: correcting errors, adding missing atoms, filling valences with hydrogens, predicting pK values for titratable amino acids, assigning predefined partial charges and radii to all atoms, and generating force field parameter/topology files for MD. Identifying, installing and effectively using the appropriate tools for each of these tasks can be difficult for novice and time-consuming for experienced users. H++ (http://biophysics.cs.vt.edu/) is a free open-source web server that automates the above key steps in the preparation of biomolecular structures for molecular modeling and simulations. H++ also performs extensive error and consistency checking, providing error/warning messages together with the suggested corrections. In addition to numerous minor improvements, the latest version of H++ includes several new capabilities and options: fix erroneous (flipped) side chain conformations for HIS, GLN and ASN, include a ligand in the input structure, process nucleic acid structures and generate a solvent box with specified number of common ions for explicit solvent MD.

Effective Born radii in the generalized Born approximation: The importance of being perfect

Generalized Born (GB) models provide, for many applications, an accurate and computationally facile estimate of the electrostatic contribution to aqueous solvation. The GB models involve two main types of approximations relative to the Poisson equation (PE) theory on which they are based. First, the self‐energy contributions of individual atoms are estimated and expressed as “effective Born radii.” Next, the atom‐pair contributions are estimated by an analytical function fGB that depends upon the effective Born radii and interatomic distance of the atom pairs. Here, the relative impacts of these approximations are investigated by calculating “perfect” effective Born radii from PE theory, and enquiring as to how well the atom‐pairwise energy terms from a GB model using these perfect radii in the standard fGB function duplicate the equivalent terms from PE theory. In tests on several biological macromolecules, the use of these perfect radii greatly increases the accuracy of the atom‐pair terms; that is, the standard form of fGB performs quite well. The remaining small error has a systematic and a random component. The latter cannot be removed without significantly increasing the complexity of the GB model, but an alternative choice of fGB can reduce the systematic part. A molecular dynamics simulation using a perfect‐radii GB model compares favorably with simulations using conventional GB, even though the radii remain fixed in the former. These results quantify, for the GB field, the importance of getting the effective Born radii right; indeed, with perfect radii, the GB model gives a very good approximation to the underlying PE theory for a variety of biomacromolecular types and conformations.

A simple clustering algorithm can be accurate enough for use in calculations of pKs in macromolecules

Structure and function of macromolecules depend critically on the ionization states of their acidic and basic groups. Most current structure‐based theoretical methods that predict pK of ionizable groups in macromolecules include, as one of the key steps, a computation of the partition sum (Boltzmann average) over all possible protonation microstates. As the number of these microstates depends exponentially on the number of ionizable groups present in the molecule, direct computation of the sum is not realistically feasible for many typical proteins that may have tens or even hundreds of ionizable groups. We have tested a simple and robust approximate algorithm for computing these partition sums for macromolecules. The method subdivides the interacting sites into independent clusters, based upon the strength of site–site electrostatic interaction. The resulting partition function is factorizable into computationally manageable components. Two variants of the approach are presented and validated on a representative test set of 602 proteins, by comparing the pK1/2 values computed by the proposed method with those obtained by the standard Monte Carlo approach used as a reference. With 95% confidence, the relative error introduced by the more accurate of the two methods is less than 0.25 pK units. The algorithms are one to two orders of magnitude faster than the Monte Carlo method, with the typical settings. A graphical representation is introduced that visualizes the clusters of strong site–site interactions in the context of the three‐dimensional (3D) structure of the macromolecule, facilitating identification of functionally important clusters of ionizable groups; the approach is exemplified on two proteins, bacteriorhodopsin and myoglobin. Proteins 2006.

Atomic level computational identification of ligand migration pathways between solvent and binding site in myoglobin

Myoglobin is a globular protein involved in oxygen storage and transport. No consensus yet exists on the atomic level mechanism by which oxygen and other small nonpolar ligands move between the myoglobins buried heme, which is the ligand binding site, and surrounding solvent. This study uses room temperature molecular dynamics simulations to provide a complete atomic level picture of ligand migration in myoglobin. Multiple trajectories—providing a cumulative total of 7 μs of simulation—are analyzed. Our simulation results are consistent with and tie together previous experimental findings. Specifically, we characterize: (i) Explicit full trajectories in which the CO ligand shuttles between the internal binding site and the solvent and (ii) pattern and structural origins of transient voids available for ligand migration. The computations are performed both in sperm whale myoglobin wild-type and in sperm whale V68F myoglobin mutant, which is experimentally known to slow ligand-binding kinetics. On the basis of these independent, but mutually consistent ligand migration and transient void computations, we find that there are two discrete dynamical pathways for ligand migration in myoglobin. Trajectory hops between these pathways are limited to two bottleneck regions. Ligand enters and exits the protein matrix in common identifiable portals on the protein surface. The pathways are located in the “softer” regions of the protein matrix and go between its helices and in its loop regions. Localized structural fluctuations are the primary physical origin of the simulated CO migration pathways inside the protein.

Analytical electrostatics for biomolecules: Beyond the generalized Born approximation

The modeling and simulation of macromolecules in solution often benefits from fast analytical approximations for the electrostatic interactions. In our previous work [G. Sigalov et al., J. Chem. Phys. 122, 094511 (2005)], we proposed a method based on an approximate analytical solution of the linearized Poisson-Boltzmann equation for a sphere. In the current work, we extend the method to biomolecules of arbitrary shape and provide computationally efficient algorithms for estimation of the parameters of the model. This approach, which we tentatively call ALPB here, is tested against the standard numerical Poisson-Boltzmann (NPB) treatment on a set of 579 representative proteins, nucleic acids, and small peptides. The tests are performed across a wide range of solvent/solute dielectrics and at biologically relevant salt concentrations. Over the range of the solvent and solute parameters tested, the systematic deviation (from the NPB reference) of solvation energies computed by ALPB is 0.5-3.5 kcal/mol, which is 5-50 times smaller than that of the conventional generalized Born approximation widely used in this context. At the same time, ALPB is equally computationally efficient. The new model is incorporated into the AMBER molecular modeling package and tested on small proteins.

Building Water Models: A Different Approach

Simplified classical water models are currently an indispensable component in practical atomistic simulations. Yet, despite several decades of intense research, these models are still far from perfect. Presented here is an alternative approach to constructing widely used point charge water models. In contrast to the conventional approach, we do not impose any geometry constraints on the model other than the symmetry. Instead, we optimize the distribution of point charges to best describe the “electrostatics” of the water molecule. The resulting “optimal” 3-charge, 4-point rigid water model (OPC) reproduces a comprehensive set of bulk properties significantly more accurately than commonly used rigid models: average error relative to experiment is 0.76%. Close agreement with experiment holds over a wide range of temperatures. The improvements in the proposed model extend beyond bulk properties: compared to common rigid models, predicted hydration free energies of small molecules using OPC are uniformly closer to experiment, with root-mean-square error <1 kcal/mol.