To request any journal reprints, please contact us.
ISIM: A Program for Grand Canonical Monte Carlo Simulations of the Ionic Environment of Biomolecules
Andreas Vitalis, Nathan A. Baker and J. Andrew McCammon
In this work we present a new software package (ISIM), which represents a flexible, computational tool for simulations of electrolyte solutions via a grand canonical Monte Carlo procedure (GCMC) with a specific capability of treating biomolecules in solution. The GCMC method provides a powerful tool for studying the ionic environments of highly charged macromolecules with attention to the atomic detail of both the solute and the mobile counterions. The ISIM software differs from previous schemes mainly by treating different ion types independently and offering a new parameterization procedure for calibrating excess chemical potentials and bulk ion concentrations. Additionally, ISIM leverages the APBS software package to provide accurate descriptions of the biomolecular electrostatic potential through the efficient solution of Poisson's equations. ISIM has been validated on a variety of test systems; we successfully reproduce elementary properties of electrolyte solutions as well as theoretical and experimental results for challenging test systems like Calmodulin and DNA.
Mechanism of Acetylcholinesterase Inhibition by Fasciculin: A 5 ns Molecular Dynamics Simulation
Kaihsu Tai, Tongye Shen, Richard H. Henchman, Yves Bourne, Pascale Marchot and J. Andrew McCammon
In our previous molecular dynamics simulation of mouse acetylcholinesterase (EC 220.127.116.11), the enzyme revealed complex gorge fluctuation. Now we report another 5 ns simulation, with acetylcholinesterase complexed with fasciculin 2. Fasciculin 2 binds to the gorge entrance of acetylcholinesterase with excellent complementarity and many polar interactions. In this simulation of the protein-protein complex,we continue to observe a two-peaked probability distribution of the gorge width, though fasciculin 2 appears to block access of ligands to the gorge. The gorge width distribution is altered when fasciculin is present: the gorge is more likely to be narrow. Though the gorge is sterically blocked by fasciculin 2, there are large increases in the opening of alternative passages, namely the side door and the back door. The catalytic triad arrangement in the acetylcholinesterase active site is disrupted with fasciculin bound. These data suggest that, in addition to the steric obstruction, fasciculin may inhibit acetylcholinesterase by allosteric and dynamical means. Additional information from these simulations can be found in reference 1 and at http://mccammon.ucsd.edu/.
Influence of Water on the Function of Acetylcholinesterase
Richard Henchman, Kaihsu Tai, Tongye Shen and J. Andrew McCammon
A 10 ns trajectory from a molecular dynamics simulation is used to examine the structure and dynamics of water around acetylcholinesterase to determine what influence water may have on its function. Particular emphasis is placed on water in the active site gorge to understand how water may affect ligand entry and binding. Despite the confining nature of the deep active site gorge ligand entry appears to be aided by a number of water properties -- fluctuations in the population of gorge waters, moderate mobility of water in the entrance to the gorge, reduced water hydrogen bonding ability, and transient cavities in the gorge. While a ligand is signifcantly slowed down by the less mobile water, this effect does not appear to be large enough to be rate limiting for enzyme catalysis of acetycholine. Images and further information are available at the website http://mccammon.ucsd.edu/.
The Displacement of Ligand through the Bottleneck Region of the Acetylcholinesterase Gorge
J. Bui and J.A. McCammon
In "Cholinesterases in the Second Millennium: Biomolecular and Pathological Aspects," N.C. Inestrosa and E.O. Campos, Eds., pp. 207-211 (2004)
Molecular dynamics simulations are applied to understand how acetylcholine can pass through the constricted region of the 20 Å deep gorge of acetylcholinesterase and gain access to the active site. The molecular dynamics (MD) simulations of the simpler ligand tetramethylammonium (TMA+) in this bottleneck were conducted using umbrella potential sampling and activated flux techniques. Potential of mean force (pmf) and barrier crossing rate constants in this constriction region were calculated. From the results of activated dynamics simulations, local conformational fluctuations of the gorge residues that were coupled to the motion of the ligand crossing were observed.
Chung F. Wong and J. Andrew McCammon
Computer-aided drug design is playing an increasing role in the development of pharmaceuticals. Computers can now be used in many ways to aid the design of therapeutic drugs. Molecular graphics software running on powerful workstations can draw useful qualitative insights from the experimental structure of a receptor and/or its complex with ligands. Rapid advances in computing technology and computational chemistry have facilitated the generation of quantitative data to help sort out promising drug leads from real and virtual chemistry libraries and to aid rational drug design. These methods range from rigid structural models employing simple but easy-to-compute scoring functions to flexible structural models with sophisticated description of intra and intermolecular interactions. Structural prediction programs have been developed to provide structural models of receptors when experimental structure is not yet available. Alternatively, useful insights into drug design can be obtained by comparative analysis of a number of inhibitors that are known to bind to a receptor. This ligand-based approach can be used without a knowledge of the three dimensional structure of the receptor. Each approach has its pros and cons and effective computer-aided drug design usually utilizes a number of methods. Here, we give a glimpse of some of these approaches with a short selected list of references that aids interested readers to further explore the subject.
Revisiting Free Energy Calculations: A Theoretical Connection to MM/PBSA and Direct Calculation of the Association Free Energy
Jessica M.J. Swanson, Richard H. Henchman and J. Andrew McCammon
The prediction of absolute ligand-receptor binding affinities is essential in a wide range of biophysical queries, from the study of protein-protein interactions to structure based drug design. End-point free energy methods, such as the Molecular Mechanics Poisson Boltzmann Surface Area (MM/PBSA) model, have received much attention and widespread application in recent literature. These methods benefit from computational efficiency as only the initial and final states of the system are evaluated, yet there remains a need for strengthening their theoretical foundation. Here a clear connection between statistical thermodynamics and end-point free energy models is presented. The importance of the association free energy, arising from one molecule's loss of translational and rotational freedom from the standard state concentration, is addressed. A novel method for calculating this quantity directly from a molecular dynamics simulation is described. The challenges of accounting for changes in the protein conformation and its fluctuations from separate simulations are discussed. A simple first order approximation of the configuration integral is presented to lay the groundwork for future efforts. This model has been applied to FKBP12, a small immunophilin that has been widely studied in the drug industry for its potential immunosuppressive and neuroregenerative effects.
The Role of Hydrogen Bonding in the Enzymatic Reaction Catalyzed by HIV-1 Protease
Joanna Trylska, Paweł Grochowski and J. Andrew McCammon
The hydrogen bond network in various stages of the enzymatic reaction catalyzed by HIV-1 protease was studied through quantum-classical molecular dynamics simulations. The approximate valence bond method was applied to the active site atoms participating directly in the rearrangement of chemical bonds. The rest of the protein with explicit solvent was treated with a classical molecular mechanics model. Two possible mechanisms were studied, general-acid/general-base (GA/GB) with Asp25 protonated at the inner oxygen, and a direct nucleophilic attack by Asp25. Strong hydrogen bonds leading to spontaneous proton transfers were observed in both reaction paths. A single-well hydrogen bond was formed between the peptide nitrogen and outer oxygen of Asp125. The proton was diffusely distributed with an average central position and transferred back and forth on a picosecond scale. In both mechanisms this interaction helped change the peptide bond hybridization, increased the partial charge on peptidyl carbon, and in the GA/GB mechanism helped deprotonate the water molecule. The inner oxygens of the aspartic dyad formed a low-barrier but asymmetric hydrogen bond; the proton was not positioned midway and made a slightly elongated covalent bond transferring from one to the other aspartate. In the GA/GB mechanism both aspartates may help deprotonate the water molecule. We observed the breakage of the peptide bond and found that the protonation of the peptidyl amine group was essential for the peptide bond cleavage. In studies of the direct nucleophilic mechanism the peptide carbon of the substrate and oxygen of Asp25 approached as close as 2.3Å.
Finite Element Solution of the Steady-state Smoluchowski Equation for Rate Constant Calculations
Yuhua Song, Yongjie Zhang, Tongye Shen, Chandrajit L. Bajaj, J. Andrew McCammon and Nathan A. Baker
This article describes the development and implementation of algorithms to study diffusion in biomolecular systems using continuum mechanics equations. Specifically, finite element methods have been developed to solve the steady-state Smoluchowski equation to calculate ligand binding rate constants for large biomolecules. The resulting software has been validated and applied to mouse acetylcholinesterase. Rates for inhibitor binding to mAChE were calculated at various ionic strengths with several different reaction criteria. The calculated rates were compared with experimental data and show very good agreement when the correct reaction criterion is used. Additionally, these finite element methods require significantly less computational resources than existing particle-based Brownian dynamics methods.
HIV-1 Protease Molecular Dynamics of a Wild-Type and of the V82F/I84V Mutant: Possible Contributions to Drug Resistance and a Potential New Target Site for Drugs
Alexander L. Perryman, Jung-Hsin Lin and J. Andrew McCammon
The protease from type 1 human immunodeficiency virus (HIV-1) is a critical drug target against which many therapeutically useful inhibitors have been developed. Because of selective pressures applied on the viral strains by the different protease inhibitors in the Highly Active Anti-Retroviral Therapies (HAART), the set of viral strains in the population has been shifting to become more drug-resistant. Because indirect effects are contributing to drug resistance, an examination of the dynamic structures of a wild type and a mutant could be insightful. Consequently, this study examined structural properties sampled during 22 nanosecond, all atom molecular dynamics (MD) simulations (in explicit water) of both a wild type and the drug-resistant V82F/I84V mutant of HIV-1 protease. The V82F/I84V mutation significantly decreases the binding affinity of all HIV-1 protease inhibitors currently used clinically. Simulations (Scott, W.R.P. & Schiffer, C.A., 2000, Structure 8, 1259-1265) have shown that the curling of the tips of the active site flaps immediately results in flap opening. In the 22ns MD simulations presented here, more frequent and more rapid curling of the mutant's active site flap tips was observed. The mutant protease's flaps opened farther than the wild type's flaps did and displayed more flexibility. This suggests that the effect of the mutations on the equilibrium between the semi-open and closed conformations could be one aspect of the mechanism of drug resistance for this mutant. In addition, correlated fluctuations in the active site and periphery were noted that point to a possible binding site for allosteric inhibitors.
Discovery of a Novel Binding Trench in HIV Integrase
Julie R. Schames, Richard H. Henchman, Jay S. Siegel, Christoph A. Sotriffer, Haihong Ni and J. Andrew McCammon
Docking of the 5CITEP inhibitor to snapshots of a 2 ns HIV-1 integrase MD trajectory indicated a previously uncharacterized trench adjacent to the active site that intermittently opens. Further docking studies of novel ligands with the potential to bind to both regions showed greater selective affinity when able to bind to the trench. Our ranking of ligands is open to experimental testing, and our approach suggests a new target for HIV-1 therapeutics.
Ribosome Motions Modulate Electrostatic Properties
Joanna Trylska, Robert Konecny, Florence Tama, Charles L. Brooks III and J. Andrew McCammon
The electrostatic properties of the 70S ribosome of Thermus thermophilus were studied qualitatively by solving the Poisson-Boltzmann (PB) equation in aqueous solution and with physiological ionic strength. The electrostatic potential was calculated for conformations of the ribosome derived by recent normal mode analysis (Tama, F. et al. Proc. Natl. Acad. Sci. USA 2003, 100, 9319--9323) of the ratchet-like reorganization that occurs during translocation (Frank, J.; Agrawal, R.K. Nature 2000, 406, 318--322). To solve the PB equation, effective parameters (charges and radii), applicable to a highly charged backbone model of the ribosome, were developed. Regions of positive potential were found at the binding site of the elongation factors G and Tu, as well as where the release factors bind. Large positive potential areas are especially pronounced around the L11 and L6 proteins. The region around the L1 protein is also positively charged, supporting the idea that L1 may interact with the E-site tRNA during its release from the ribosome after translocation. Functional rearrangement of the ribosome leads to electrostatic changes which may help the translocation of the tRNAs during the elongation stage.
PDB2PQR: An Automated Pipeline for the Setup, Execution and Analysis of Poisson-Boltzmann Electrostatics Calculations
Todd J. Dolinsky, Jens E. Nielsen, J. Andrew McCammon and Nathan A. Baker
Continuum solvation models, such as Poisson-Boltzmann and Generalized Born methods, have become increasingly popular tools for investigating the influence of electrostatics on biomolecular structure, energetics and dynamics. However, the use of such methods requires accurate and complete structural data as well as force field parameters such as atomic charges and radii. Unfortunately, the limiting step in continuum electrostatics calculations is often addition of missing atomic coordinates to molecular structures from the Protein Data Bank and the assignment of parameters to biomolecular structures. To address this problem, we have developed the PDB2PQR web service (http://www.poissonboltzmann.org/pdb2pqr). This server automates many of the common tasks of preparing structures for continuum electrostatics calculations, including: adding a limited number of missing heavy atoms to biomolecular structures, estimating titration states and protonating biomolecules in a manner consistent with favorable hydrogen-bonding, assigning charge and radius parameters from a variety of force fields, and finally generating "PQR" output compatible with several popular computational biology packages. This service is intended to facilitate the setup and execution of electrostatics calculations for both experts and non-experts and thereby broaden the accessibility of the biological community to continuum electrostatics analyses of biomolecular systems.
Accelerated Molecular Dynamics: A Promising and Efficient Simulation Method for Biomolecules
Donald Hamelberg, John Mongan and J. Andrew McCammon
Many interesting dynamic properties of biological molecules cannot be simulated directly using molecular dynamics because of nanosecond timescale limitations. These systems are trapped in potential energy minima with high free energy barriers for large numbers of computational steps. The dynamic evolution of many molecular systems occurs through a series of rare events as the system moves from one potential energy basin to another. Therefore, we have proposed a robust bias potential function that can be used in an efficient accelerated molecular dynamics approach to simulate the transition of high energy barriers without any advance knowledge of the location of either the potential energy wells or saddle points. In this method, the potential energy landscape is altered by adding a bias potential to the true potential such that the escape rates from potential wells are enhanced, which accelerates and extends the time scale in molecular dynamics simulations. Our definition of the bias potential echoes the underlying shape of the potential energy landscape on the modified surface, thus allowing for the potential energy minima to be well defined, and hence properly sampled during the simulation. We have shown that our approach, which can be extended to biomolecules, samples the conformational space more efficiently than normal molecular dynamics simulations, and converges to the correct canonical distribution.
Acetylcholinesterase: Enhanced Fluctuations and Alternative Routes to the Active Site in the Complex with Fasciculin-2
Jennifer M. Bui, Kaihsu Tai and J. Andrew McCammon
A 15ns molecular dynamics simulation is reported for the complex of mouse acetylcholinesterase (mAChE) and the protein neurotoxin fasciculin-2. Compared to a 15ns simulation of apo-mAChE, the structural fluctuations of the enzyme are substantially increased in magnitude for the enzyme in the complex. Fluctuations of part of the long omega loop (residues 69-96) are particularly enhanced. This loop forms one wall of the active site, and the enhanced fluctuations lead to additional routes of access to the active site.
Charge Optimization of the Interface between Protein Kinases and their Ligands
Peter A. Sims, Chung F. Wong and J. Andrew McCammon
Examining the potential for electrostatic complementarity between a ligand and a receptor is a useful technique for rational drug design and can demonstrate how a system prioritizes interactions when allowed to optimize its charge distribution. In this computational study, we implemented the previously developed, continuum solvent-based charge optimization theory with a simple, quadratic programming algorithm and the UHBD Poisson-Boltzmann solver. This method allows one to compute the best set of point charges for a ligand or ligand region based on the ligand and receptor shape, and the receptor partial charges, by optimizing the binding free energy obtained from a continuum-solvent model. We applied charge optimization to a fragment of the heat-stable protein kinase inhibitor (PKI) of protein kinase A (PKA), to three flavopiridol inhibitors of CDK2, and to cyclin A which interacts with CDK2 to regulate the cell cycle. We found that a combination of global (involving every charge) and local (involving only charges in a local region) optimization can give useful hints for designing better inhibitors. Although some parts of an inhibitor may already contribute significantly to binding, we found that they could still be the most important targets for modifications to obtain stronger binders. In studying the binding of flavopiridol inhibitors to CDK2, comparable binding affinity could be obtained regardless of whether the net charges of the inhibitors were constrained to -2, -1, 0, 1, or 2 during the optimization. This provides flexibility in inhibitor design when a certain net charge of the inhibitor is desired in addition to strong binding affinity. For the study of the PKA-PKI and CDK2-cylin A interfaces, we identified residues whose charge distributions are already close to optimal and those whose charge distributions could be refined to further improve binding.
Standard Free Energy of Releasing a Localized Water Molecule from the Binding Pockets of Proteins: Double-Decoupling Method
Donald Hamelberg and J. Andrew McCammon
Localized water molecules in the binding pockets of proteins play an important role in non-covalent association of proteins and small drug compounds. At times the dominant contribution to the binding free energy comes from the release of localized water molecules in the binding pockets of biomolecules. Therefore, in order to quantify the energetic importance of these water molecules for drug design purposes, we have used the double-decoupling approach to calculate the standard free energy of tying up a water molecule in the binding pockets of two protein complexes. The double-decoupling approach is based on the underlying principle of statistical thermodynamics. We have calculated the standard free energies of tying up the water molecule in the binding pockets of these complexes to be favorable. These water molecules stabilize the protein-drug complexes by interacting with the ligands and binding pockets of the protein complexes. Our results offer ideas that could be used in optimizing protein-drug interactions, by designing ligands that are capable of targeting localized water molecules in the binding pockets of proteins. The resulting free energy of ligand binding could benefit from the potential free energy gain accompanying the release of these water molecules. Furthermore, we have examined the theoretical background of the double-decoupling method and its connection to the molecular dynamics thermodynamic integration techniques.
Electrostatic Interaction between RNA and Protein Capsid in CCMV Simulated by a Coarse-grain RNA model and a Monte Carlo Approach
Deqiang Zhang, Robert Konecny, Nathan A. Baker and J. Andrew McCammon
Although many viruses have been crystallized and the protein capsid structures have been determined by x-ray crystallography, the nucleic acids often cannot be resolved. This is especially true for RNA viruses. The lack of information about the conformation of DNA/RNA greatly hinders our understanding of the assembly mechanism of various viruses. Here we combine a coarse-grain model and a Monte Carlo method to simulate the distribution of viral RNA inside the capsid of cowpea chlorotic mottle virus. Our results show that there is very strong interaction between the N-terminal residues of the capsid proteins, which are highly positive charged, and the viral RNA. Without these residues, the binding energy disfavors the binding of RNA by the capsid. The RNA forms a shell close to the capsid with the highest densities associated with the capsid dimers. These high-density regions are connected to each other in the shape of a continuous net of triangles. The overall icosahedral shape of the net overlaps with the capsid subunit icosahedral organization. Medium density of RNA is found under the pentamers of the capsid. These findings are consistent with experimental observations.
Finite concentration effects on diffusion-controlled reactions
Sanjib Senapati, Chung F. Wong and J. Andrew McCammon
The algorithm by Northrup, Allison, and McCammon [J. Chem. Phys. 80, 1517 (1984)] has been used for two decades for calculating the diffusion-influenced rate-constants of enzymatic reactions. Although many interesting results have been obtained, the algorithm is based on the assumption that substrate/substrate interactions can be neglected. This approximation may not be valid when the concentration of the ligand is high. In this work, we constructed a simulation model that can take substrate/substrate interactions into account. We first validated the model by carrying out simulations in ways that could be compared to analytical theories. We then carried out simulations to examine the possible effects of substrate/substrate interactions on diffusion-controlled reaction rates. For a substrate concentration of 0.1 mM, we found that the diffusion-controlled reaction rates were not sensitive to whether substrate/substrate interactions were included. On the other hand, we observed significant influence of substrate/substrate interactions on calculated reaction rates at a substrate concentation of 0.1M. Therefore, a simulation model that takes substrate/substrate interactions into account is essential for reliably predicting diffusion-controlled reaction rates at high substrate concentrations, and one such simulation model is presented here.
Constant pH Molecular Dynamics in Generalized Born Implicit Solvent
John Mongan, David A. Case and J. Andrew McCammon
A new method is proposed for constant pH molecular dynamics (MD), employing generalized Born (GB) electrostatics. Protonation states are modeled with different charge sets, and titrating residues sample a Boltzmann distribution of protonation states as the simulation progresses, using Monte Carlo sampling based on GB-derived energies. The method is applied to four different crystal structures of hen egg-white lysozyme (HEWL). pKa predictions derived from the simulations have root-mean-square (RMS) error of 0.82 relative to experimental values. Similarity of results between the four crystal structures shows the method to be independent of starting crystal structure; this is in contrast to most electrostatics-only models. A strong correlation between conformation and protonation state is noted and quantitatively analyzed, emphasizing the importance of sampling protonation states in conjunction with dynamics.