To request any journal reprints, please contact us.
Small Molecule pKa Prediction with Continuum Electrostatics Calculations
Michael J. Potter, Michael K. Gilson and J. Andrew McCammon
In the simulation of biomolecules and their interactions, it is critically important to correctly model the protonation states of titratable groups. Ideally, such models would allow for adjustment of the protonation states to reflect conformational and other environmental changes. It is thus desirable to develop theoretical methods capable of quickly and accurately predicting pKa values. A number of technique based on Poisson-Boltzmann electrostatics calculations have recently been described for application to proteins. In the method of Antosiewicz et al., each group is assigned an "initial" pKa which represents its behavior when the group is isolated in solution. Finite-difference solutions to the linearized Poisson-Boltzmann equation are then used to adjust this pKa to reflect the group's electrostatic environment within the protein. To test this approach further and to explore its utility for other types of molecules, we have applied the method to a variety of experimentally well-characterized diamines and diacids. The compounds considered were ethylenediamine, 1,3-diaminopropane, 1,2-diaminopropane, 1,4-diaminobutane, malonic acid, succinic acid, and glutaric acid.
Biological Applications of Electrostatic Calculations and Brownian Dynamics Simulations
Jeffry D. Madura, Malcolm E. Davis, Michael K. Gilson, Rebecca C. Wade, Brock A. Luty and J. Andrew McCammon
Extended Hellmann-Feynman Theorem for Non-Stationary States and Its Application in Quantum-Classical Molecular Dynamics Simulations
P. Bała, B. Lesyng and J.A. McCammon
An extended Hellmann-Feynman theorem for non-stationary states is presented. For systems in which some particles are described quantum-mechanically and others classically, the extended theorem allows one to determine more precisely the forces exerted by quantum particles on classical particles. The forces are used in quantum-classical molecular dynamics (QCMD) simulations, where the extended theorem allows an increase in the integration time-step and a reduction in the energy drift. The QCMD model also accounts for a dynamical vibronic coupling between quantum and classical subsystems, which further improves the quality of the simulation results.
Simulation of Enzyme-Substrate Encounter with Gated Active Sites
Rebecca C. Wade, Brock A. Luty, Eugene Demchuk, Jeffry D. Madura, Malcolm E. Davis, James M. Briggs and J. Andrew McCammon
We describe a brownian dynamics simulation method that allows investigation of the effects of receptor flexibility on ligand binding rates. The method is applied to the encounter of substrate, glyceraldehyde 3-phosphate, with triose phosphate isomerase, a diffusion-controlled enzyme with flexible peptide loops at its active sites. The simulations show that while the electrostatic field surrounding the enzyme steers the substrate into its active sites, the flexible loops appear to have little influence on the substrate binding rate. The dynamics of the loops may therefore have been optimized during evolution to minimize their interference with the substrate's access to the active sites. The calculated and experimental rate constants are in good agreement.
Applications of Quantum-Classical and Quantum-Stochastic Molecular Dynamics Simulations for Proton Transfer Processes
P. Bała, B. Lesyng and J.A. McCammon
Quantum-classical and quantum-stochastic molecular dynamics models (QCMD/QSMD) are formulated and applied to describe proton transfer processes in three model systems - the proton bound ammonia-ammonia dimer in an external electrostatic field; malonaldehyde, which undergoes a quantum tautomeric rearrangement; and phospholipase A2, an enzyme which induces a water dissociation process in its active site followed by proton hopping to a histidine imidazole ring. The proton dynamics are described by the time-dependent Schrödinger equation. The dynamics of the classical atoms are described using classical molecular dynamics. Coupling between the quantum proton (s) and the classical atoms is accomplished via conventional or extended Hellmann-Feynman forces, as well as the time-dependence of the potential energy function in the Schrödinger equation. The interaction of the system with its environment is described by stochastic forces. Possible extensions of the models as well as future applications in molecular structure and dynamics analysis will be briefly discussed.
Ordinary Differential Equations of Molecular Dynamics
J.A. McCammon, B.M. Pettitt and L.R. Scott
The ordinary differential equations of Newtonian dynamics are used in atomic simulations with the method of molecular dynamics. The basic issues are surveyed and standard algorithms are described. Several algorithmic variants are discussed. Some advanced ideas relating to parallel computation are considered.
Separation-shifted Scaling, a new Scaling Method for Lennard-Jones Interactions in Thermodynamic Integration
M. Zacharias, T.P. Straatsma and J.A. McCammon
A new method of simultaneously scaling and shifting the Lennard-Jones (LJ) potential in molecular dynamics (MD) and thermodynamic integration (TI) simulations is presented. The approach allows the smooth creation or annihilation of atoms or molecules in an ensemble of solvent molecules during a molecular simulation. By scaling and shifting the LJ potential in the direction of the interatomic distance between particles, the method eliminates the problem of the creation or annihilation of a large repulsive LJ potential at the initial or final state of a TI. The optimal degree of shifting and scaling the LJ potential as a function of a control variable λ was studied for the annihilation and creation of neon in aqueous solution. The method was further tested on the calculation of the free energy of aqueous solvation of a small molecule, ethanol. In contrast to linear scaling of the LJ potential during TI, the calculated free energies using the new separation-shifted scaling approach are reasonably well converged after 200-500 ps of simulation and show smaller hysteresis comparing forward and reverse TI.
Molecular Modeling Methods. Basic Techniques and Challenging Problems
Bogdan Lesyng and J. Andrew McCammon
An overview is presented of computer modeling and simulation methods that play an increasing role in drug design: quantum chemical methods, molecular mechanics, molecular dynamics and Brownian dynamics. The application of molecular dynamics for the prediction of thermodynamic properties like free energy differences and binding constants is discussed. The Brownian dynamics method is presented in connection with the calculation of effective electrostatic forces using the Poisson-Boltzmann equation, which allows one to sample ligand-binding geometries and to predict the kinetics of diffusion-limited enzyme reactions. New techniques that have recently been extensively developed, such as the global energy minimization and quantum-classical dynamics methods, are also introduced. The molecular modeling methods are illustrated with selected examples.
Open "Back Door" in a Molecular Dynamics Simulation of Acetylcholinesterase
M.K. Gilson, T.P. Straatsma, J.A. McCammon, D.R. Ripoll, C.H. Faerman, P.H. Axelsen, I. Silman and J.L. Sussman
The enzyme acetylcholinesterase generates a strong electrostatic field that can attract the cationic substrate acetylcholine to the active site. However, the long and narrow active site gorge seems inconsistent with the enzyme's high catalytic rate. A molecular dynamics simulation of acetylcholinesterase in water reveals the transient opening of a short channel, large enough to pass a water molecule, through a thin wall of the active site near tryptophan-84. This simulation suggests that substrate, products, or solvent could move through this "back door," in addition to the entrance revealed by the crystallographic structure. Electrostatic calculations show a strong field at the back door, oriented to attract the substrate and the reaction product choline and to repel the other reaction product, acetate. Analysis of the open back door conformation suggests a mutation that could seal the back door and thus test the hypothesis that thermal motion of this enzyme may open multiple routes of access to its active site.
Prediction of pH-Dependent Properties of Proteins
Jan Antosiewicz, J. Andrew McCammon and Michael K. Gilson
We describe what may be the most accurate approach currently available for the calculation of the pKas of ionizable groups in proteins. The accuracy is assessed by comparison of computed pKas with 60 measured pKas in a total of seven proteins. The overall root-mean-square error is 0.89 pKa units. Linear regression analysis of computed versus measured pKas yields a slope of 0.95, y-intercept of -0.02 and a correlation coefficient of 0.96. The proposed approach also picks out many of the shifted pKas of groups in enzyme active sites and special salt bridges. However, it does yield several over-shifted pKas and tends to underestimate pKa shifts which result from desolvation effects. We examine the ability of the new approach to reproduce the dependence of protein stability upon Ph, using the ionization polynomial formalism. Overall features of the stability curves are reproduced, but the quantitative agreement is not particularly good. The reasons for the disagreement may have to do both with insufficient accuracy in the theory and with uncertainty in the nature of the unfolded state of proteins. The methodology described here is based upon finite difference solutions of the Poisson-Boltzmann equation. Its success depend upon the use of the rather high protein dielectric constant of 20. However, theoretical considerations and the fact that pKa shifts which result from desolvation are underestimated here imply that the dielectric constant of the protein interior actually is lower than 20. We suggest that the high protein dielectric constant improves the overall agreement with experiment because it accounts approximately for phenomena which tend to mitigate pKa shifts and which are not specifically included in the model. These include conformational relaxation and specific ion-binding. Future models based upon a low protein dielectric constant and treating such phenomena explicitly might yield improved agreement with experiment.
Combined Conformational Search and Finite-Difference Poisson-Boltzmann Approach for Flexible Docking. Application to an Operator Mutation in the Lambda Repressor-Operator Complex
M. Zacharias, B.A. Luty, M.E. Davis and J.A. McCammon
The N-terminal domain of the phage λ repressor binds as a dimer to its palindromic DNA operator sequence. In addition to a helix-turn-helix DNA recognition motif, the first six amino acids of the phage λ repressor form a flexible peptide segment which wraps around DNA. Site-directed mutagenesis studies have shown that amino acid replacements or partial removal of the arm structure, or changes in the DNA sequence contacting the N-terminal arm, can lower the repressor-operator binding affinity by several orders of magnitude. The finite-difference Poisson-Boltzmann approach in combination with a conformational search procedure was used to study energetic contributions of the λ arm to repressor-operator recognition based on the high resolution X-ray structure. It allows for the local relaxation of the structure upon changing the DNA sequence in the λ arm binding region. A simplified potential energy function including torsional, truncated Lennard-Jones and approximate electrostatic terms is used in the initial step to screen out energetically unfavorable structures. The electrostatic energy of selected conformations is subsequently calculated more accurately using the finite-difference Poisson-Boltzmann approach. The method was applied to study the effect of a C→T mutation at position 6 of the consensus half-site of the operator. This base-pairs contacts Lys4 which is part of the arm segment. Keeping only the Lys4 side-chain mobile and with the wild-type DNA operator sequence, several conformations close to the X-ray structure were identified as those with lowest energy. In the case of the DNA mutation, lowest energy conformations differed significantly from those selected for the wild-type sequence. These initial calculations indicate that the approach might be a useful tool to estimate conformational and energetic effects upon mutagenesis of protein-DNA complexes.
Acetylcholinesterase: Effects of Ionic Strength and Dimerization on the Rate Constants
J. Antosiewicz, M.K. Gilson and J.A. McCammon
Israel Journal of Chemistry, Vol. 34, pp. 151-158 (1994)
Free Energy and Binding Selectivity
J. Andrew McCammon
A brief summary is presented of a discussion on the use of free energy simulations to interpret and predict molecular recognition.
Sulfate Anion in Water: Model Structural, Thermodynamic, and Dynamic Properties
William R. Cannon, B. Montgomery Pettitt and J. Andrew McCammon
Binding energies, intermolecular distances, and partial charges from ab initio studies of SO42- + H2O were used to develop a molecular mechanics model for SO42- in water. Structural, dynamic, and thermodynamic properties for the resulting model used in condensed-phase simulations agree well with experimental estimates. Thirteen waters are found to be present in the first solvation shell, and the residence time of these waters has been calculated to be 23 ps. The model anion has a free energy of solvation relative to xenon of -275 kcal/mol, which compares well with the experimental estimate of -266 kcal/mol.
Parallelization using Spatial Decomposition for Molecular Dynamics
Terry W. Clark, Reinhard von Hanxleden, J. Andrew McCammon and L. Ridgway Scott
Several algorithms have been used for parallel molecular dynamics, including the replicated algorithm and those based on spatial decompositions. The replicated algorithm stores the entire system's coordinates and forces at each processor, and therefore has a low overhead in maintaining the data distribution. Spatial decompositions distribute the data, providing better locality and scalability with respect to memory and computation. We present EULERGROMOS, a parallelization of the GROMOS molecular dynamics program which is based on a spatial decomposition. EULERGROMOS parallelizes all molecular dynamics phases, with most data structures using O(N/P) memory. The paper focuses on the structure of EULERGROMOS and analyses its performance using molecular systems of current interest in the molecular dynamics community. EULERGROMOS achieves performance increases with as few as twenty atoms per processor. We also compare EULERGROMOS with an earlier parallelization of GROMOS, UHGROMOS, which uses the replicated algorithm.
Pepzyme Dynamics and Conformation: A Molecular Dynamics Study in Water
Tami J. Marrone and J. Andrew McCammon
Recently, Atassi and Manshouri described two 29-residue cyclic peptides that were designed to perform the same function as the enzymes trypsin and chymotrypsin. Residues of the corresponding active site were linked using glycine spacers to maintain certain distances between residues. The cyclic peptide was closed by making a disulfide bridge between the N- and C-terminus cysteines. The trypsin-based "pepzyme," TrPepz, is reported to have the same specificity as trypsin but a somewhat lower catalytic rate constant. We have examined TrPepz through molecular dynamics simulations to provide some insight into the flexibility and chemistry of this molecule. Simulations in water of the free pepzyme and its complex with lysine show that both systems are quite flexible, especially the free pepzyme. The reported binding and catalytic properties cannot be reconciled with the simulation results unless, perhaps, one assumes that a number of distinct conformations are functionally competent.
Treatment of Rotational Isomeric States III. The Use of Biasing Potentials
T.P. Straatsma and J.A. McCammon
A technique is described to determine umbrella biasing potentials that can be used to enhance sampling of rotational isomeric states in molecular simulations of polypeptides in water. The analytical biasing potential functions are obtained through fitting of potentials of mean force obtained by thermodynamic integration simulations to a small number of functions used to describe dihedral torsion. The resulting dramatic increase in efficiency of sampling is illustrated by comparison of molecular simulations of the alanine-dipeptide molecule in aqueous solution, with and without the use of the biasing potentials. The same biasing potentials were used in simulations of the alanine-tripeptide and the alanine-heptapeptide in aqueous solution. Similar dramatic increases in sampling efficiency were observed for these simulations, when the biasing potentials were applied, which suggests that these biasing potentials, although determined for the dipeptide, may be transferable to larger peptide chains. It is illustrated how the biasing potentials can be used in biasing potential annealing simulations, leading to rapid folding of peptide chains into aqueous solution structures with low energy.