New-generation amber united-atom force field.

We have developed a new-generation Amber united-atom force field for simulations involving highly demanding conformational sampling such as protein folding and protein-protein binding. In the new united-atom force field, all hydrogens on aliphatic carbons in all amino acids are united with carbons except those on Calpha. Our choice of explicit representation of all protein backbone atoms aims at minimizing perturbation to protein backbone conformational distributions and to simplify development of backbone torsion terms. Tests with dipeptides and solvated proteins show that our goal is achieved quite successfully. The new united-atom force field uses the same new RESP charging scheme based on B3LYP/cc-pVTZ//HF/6-31g** quantum mechanical calculations in the PCM continuum solvent as that in the Duan et al. force field. van der Waals parameters are empirically refitted starting from published values with respect to experimental solvation free energies of amino acid side-chain analogues. The suitability of mixing new point charges and van der Waals parameters with existing Amber covalent terms is tested on alanine dipeptide and is found to be reasonable. Parameters for all new torsion terms are refitted based on the new point charges and the van der Waals parameters. Molecular dynamics simulations of three small globular proteins in the explicit TIP3P solvent are performed to test the overall stability and accuracy of the new united-atom force field. Good agreements between the united-atom force field and the Duan et al. all-atom force field for both backbone and side-chain conformations are observed. In addition, the per-step efficiency of the new united-atom force field is demonstrated for simulations in the implicit generalized Born solvent. A speedup around two is observed over the Duan et al. all-atom force field for the three tested small proteins. Finally, the efficiency gain of the new united-atom force field in conformational sampling is further demonstrated with a well-known toy protein folding system, an 18 residue polyalanine in distance-dependent dielectric. The new united-atom force field is at least a factor of 200 more efficient than the Duan et al. all-atom force field for ab initio folding of the tested peptide.

Chelate effect in cyclodextrin dimers: a computational (MD, MM/PBSA, and MM/GBSA) study.

The complexation of an adamantyl-phosphate derivative with one beta-cyclodextrin, with two beta-cyclodextrins, and with two beta-cyclodextrins dimerized with a disulfide bridge was studied by computational methods (MD, MM/PBSA, and MM/GBSA) to analyze and rationalize the chelate effect. Although this effect is usually explained by invoking favorable entropy contribution due to the preorganization of the ligand, it has been determined experimentally that in this case it is enthalpy-driven. The computational results are in accord with this finding, although the entropy contribution due to the solvent structural organization around the complex is crucial for the final estimates of the free energy of complexation.

Automatic atom type and bond type perception in molecular mechanical calculations.

In molecular mechanics (MM) studies, atom types and/or bond types of molecules are needed to determine prior to energy calculations. We present here an automatic algorithm of perceiving atom types that are defined in a description table, and an automatic algorithm of assigning bond types just based on atomic connectivity. The algorithms have been implemented in a new module of the AMBER packages. This auxiliary module, antechamber (roughly meaning "before AMBER"), can be applied to generate necessary inputs of leap-the AMBER program to generate topologies for minimization, molecular dynamics, etc., for most organic molecules. The algorithms behind the manipulations may be useful for other molecular mechanical packages as well as applications that need to designate atom types and bond types.

Hierarchical database screenings for HIV-1 reverse transcriptase using a pharmacophore model, rigid docking, solvation docking, and MM-PB/SA.

In this work, an efficient strategy was presented to search drug leads for human immunodeficiency virus type 1 reverse transcriptase (HIV-1 RT) using hierarchical database screenings, which included a pharmacophore model, multiple-conformation rigid docking, solvation docking, and molecular mechanics-Poisson-Boltzmann/surface area (MM-PB/SA) sequentially. Encouraging results were achieved in searching a refined available chemical directory (ACD) database: the enrichment factor after the first three filters was estimated to be 25-fold; the hit rate for all the four filters was predicted to be 41% in a control test using 37 known HIV-1 non-nucleoside reverse transcriptase inhibitors; 10 out of 30 promising solvation-docking hits had MM-PB/SA binding free energies better than -6.8 kcal/mol and the best one, HIT15, had -17.0 kcal/mol. In conclusion, the hierarchical multiple-filter database searching strategy is an attractive strategy in drug lead exploration.

Prediction of pKa shifts in proteins using a combination of molecular mechanical and continuum solvent calculations.

The prediction of pKa shifts of ionizable groups in proteins is of great relevance for a number of important biological phenomena. We present an implementation of the MM-GBSA approach, which combines molecular mechanical (MM) and generalized Born (GB) continuum solvent energy terms, to the calculation of pKa values of a panel of nine proteins, including 69 individual comparisons with experiment. While applied so far mainly to the calculation of biomolecular binding free energies, we show that this method can also be used for the estimation of protein pKa shifts, with an accuracy around 1 pKa unit, even for strongly shifted residues. Our analysis reveals that the nonelectrostatic terms that are part of the MM-GBSA free energy expression are important contributors to improved prediction accuracy. This suggests that most of the previous approaches that focus only on electrostatic interactions could be improved by adding other nonpolar energy terms to their free energy expression. Interestingly, our method yields best accuracy at protein dielectric constants of epsilonint = 2-4, which is in contrast to previous approaches that peak at higher epsilonint > or = 8. An important component of our procedure is an intermediate minimization step of each protonation state involving different rotamers and tautomers as a way to explicitly model protein relaxation upon (de)protonation.

Development and testing of a general amber force field.

We describe here a general Amber force field (GAFF) for organic molecules. GAFF is designed to be compatible with existing Amber force fields for proteins and nucleic acids, and has parameters for most organic and pharmaceutical molecules that are composed of H, C, N, O, S, P, and halogens. It uses a simple functional form and a limited number of atom types, but incorporates both empirical and heuristic models to estimate force constants and partial atomic charges. The performance of GAFF in test cases is encouraging. In test I, 74 crystallographic structures were compared to GAFF minimized structures, with a root-mean-square displacement of 0.26 A, which is comparable to that of the Tripos 5.2 force field (0.25 A) and better than those of MMFF 94 and CHARMm (0.47 and 0.44 A, respectively). In test II, gas phase minimizations were performed on 22 nucleic acid base pairs, and the minimized structures and intermolecular energies were compared to MP2/6-31G* results. The RMS of displacements and relative energies were 0.25 A and 1.2 kcal/mol, respectively. These data are comparable to results from Parm99/RESP (0.16 A and 1.18 kcal/mol, respectively), which were parameterized to these base pairs. Test III looked at the relative energies of 71 conformational pairs that were used in development of the Parm99 force field. The RMS error in relative energies (compared to experiment) is about 0.5 kcal/mol. GAFF can be applied to wide range of molecules in an automatic fashion, making it suitable for rational drug design and database searching.

Enhanced ab initio protein folding simulations in Poisson-Boltzmann molecular dynamics with self-guiding forces.

We have investigated the sampling efficiency in molecular dynamics with the PB implicit solvent when self-guiding forces are added. Compared with a high-temperature dynamics simulation, the use of self-guiding forces in room-temperature dynamics is found to be rather efficient as measured by potential energy fluctuation, gyration radius fluctuation, backbone RMSD fluctuation, number of unique clusters, and distribution of low RMSD structures over simulation time. Based on the enhanced sampling method, we have performed ab initio folding simulations of two small proteins, betabetaalpha1 and villin headpiece. The preliminary data for the folding simulations is presented. It is found that betabetaalpha1 folding proceeds by initiation of the turn and the helix. The hydrophobic collapse seems to be lagging behind or at most concurrent with the formation of the helix. The hairpin stability is weaker than the helix in our simulations. Its role in the early folding events seems to be less important than the more stable helix. In contrast, villin headpiece folding proceeds first by hydrophobic collapse. The formation of helices is later than the collapse phase, different from the betabetaalpha1 folding.

Investigation of neuraminidase-substrate recognition using molecular dynamics and free energy calculations.

Development of the new generation of therapeutics against the influenza viral coat protein neuraminidase is a response to the continuing threat of influenza epidemics. A variety of structurally similar compounds have been reported that vary greatly in their ability to inhibit neuraminidase, a critical enzyme that cleaves sialic acid and promotes virion release. To determine how neuraminidase exhibits this wide range of affinities with structurally similar compounds, molecular dynamic simulations, coupled with free energy calculations, were used to determine the binding components of a series of neuraminidase inhibitors. Using four cocrystal structures of neuraminidase-inhibitor complexes, we examined the structural and energetic components of ligand potency and selectivity. An in-depth energetic analysis, including internal energy, entropy, and nonbonded interactions, reveals that potency of ligand binding is governed by nonpolar contacts. Electrostatic components generally oppose binding, although two of the best inhibitors use electrostatic interactions to orient the ligand. This investigation suggests that the enhanced selectivity and potency of the better ligands may arise from an improved positioning of their ligand atoms in the active site due to polar and hydrophobic functionalities. Simulations that included crystal water molecules in the active site indicate that the more potent ligands make less use of water-mediated interactions.

Direct hydroxide attack is a plausible mechanism for amidase antibody 43C9.

Direct hydroxide attack on the scissile carbonyl of the substrate has been suggested as a likely mechanism for esterase antibodies elicited by phosphonate haptens, which mimic the transition states for the alkaline hydrolysis of esters.1 The unique amidase activity of esterase antibody 43C9 has been attributed to nucleophilic attack by an active-site histidine residue.2 Yet, the active site of 43C9 is strikingly similar to those of other esterase antibodies, particularly 17E8. We have carried out quantum mechanical calculations, molecular dynamics simulations, and free energy calculations to assess the mechanism involving direct hydroxide attack for 43C9. Results support this mechanism and suggest that the mechanism is plausible for other antiphosphonate antibodies that catalyze the hydrolysis of (p-nitro)phenyl esters.

Closing of the nucleotide pocket of kinesin-family motors upon binding to microtubules.

We have used adenosine diphosphate analogs containing electron paramagnetic resonance (EPR) spin moieties and EPR spectroscopy to show that the nucleotide-binding site of kinesin-family motors closes when the motor.diphosphate complex binds to microtubules. Structural analyses demonstrate that a domain movement in the switch 1 region at the nucleotide site, homologous to domain movements in the switch 1 region in the G proteins [heterotrimeric guanine nucleotide-binding proteins], explains the EPR data. The switch movement primes the motor both for the free energy-yielding nucleotide hydrolysis reaction and for subsequent conformational changes that are crucial for the generation of force and directed motion along the microtubule.

Free energy calculations for theophylline binding to an RNA aptamer: Comparison of MM-PBSA and thermodynamic integration methods.

We have applied the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method (J. Srinivasan, T. E. Cheatham, P. Cieplak, P. A. Kollman, and D. A. Case, Journal of the American Chemical Society, 1998, Vol. 120, pp. 9401-9409) to study the interaction of an RNA aptamer with theophylline and its analogs. The MM-PBSA free energy analysis provides a reasonable absolute binding free energy for the RNA aptamer-theophylline complex formation. Energetic analysis reveals that the van der Waals interaction and the nonpolar contribution to solvation provide the basis for the favorable absolute free energy of complex. This trend is similar to other protein-ligand interactions studied previously. The MM-PBSA method also ranks the relative binding energies of five theophylline analogs approximately correctly, but not as well as the more conventional thermodynamic integration calculations, which were carried out to convert theophylline into its analogs. The comparison of MM-PBSA with TI suggests that the MM-PBSA method has some difficulties with the first-solvation-shell energetics.

The intramolecular mechanism for the second proton transfer in triosephosphate isomerase (TIM): a QM/FE approach.

The intramolecular mechanism we earlier proposed [Alagona, G.; Desmeules, P.; Ghio, C.; Kollman, P. A. J Am Chem Soc 1984, 106, 3623] for the second proton transfer of the reaction catalyzed by triosephosphate isomerase (TIM) is examined ab initio at the HF and MP2/6-31+G** levels in vacuo for two conformers of the enediolate phosphate (ENEP), with the ethereal oxygen of the phosphate group either syn (X), as in the crystal structure, or anti (Y) with respect to the enediolate carbonyl O. The barrier height for the intramolecular proton transfer occurring in enediolate is very sensitive to electron correlation corrections. The MP2 internal energy barrier is much lower than the HF one, while the free energy (FE) barrier is even more favorable, indicating that the enzyme presence is not requested to speed up that step. An investigation of the dynamical aspects of the mechanism, along the pathway from ENEP A (with H on O(1)) to TS and from TS to ENEP B (with H on O(2)), was, however, carried out in the presence of the enzyme field while using a neutral His-95 with its proton on Ndelta. To perform the FE simulations, it was necessary to parametrize in the AMBER force-field the ENEP A, TS and B species, whose partial charges have been determined with the RESP procedure, with the X and Y arrangements of the phosphate head. Actually, the FE/QM approach produced a low barrier and a substantial balance between A and B, especially at the MP2 level. The trajectories, analyzed paying a particular attention to the positions assumed by His-95 and by the other active site residues, put forward somewhat different H-bond patterns around the X or Y enediolate phosphate.

Molecular dynamics simulation study of the negative correlation in antibody AZ28-catalyzed oxy-cope rearrangement.

The oxy-Cope rearrangement reaction in the antibody AZ28 is investigated using ab initio molecular orbital calculations and molecular mechanical molecular dynamics simulations. This antibody, AZ28, is known as one of the few systems where the mature catalytic antibody shows a negative correlation between the transition state analogue (TSA) binding affinity and the catalytic rate of the oxy-Cope rearrangement compared to the germ line catalytic antibody. The ab initio optimized structure shows that the transition state structure has a more planar configuration than the TSA. The favorable electrostatic interactions between AZ28 and the transition state analogue overcome the unfavorable van der Waals interactions; thus, AZ28 shows higher binding affinity for the TSA than the germ line. However, the AZ28 is not flexible enough to accept the relatively planar transition state structure. Because the lower flexibility causes poorer antibody-hapten interaction energies, the activation free energy of the oxy-Cope rearrangement becomes larger in the mature antibody than the germ line. We show that the differences in flexibility between the germ line and the mature form and the differences in structure between TSA and the transition state are the origin of the negative correlation in AZ28-catalyzed oxy-Cope rearrangement. The mutation of residue 34 of the light chain, 34(L), affects the binding free energies through the interresidue interaction because it is the closest to the hapten among the six mutatable residues. However, it does not affect the negative correlation.

pKa, MM, and QM studies of mechanisms of beta-lactamases and penicillin-binding proteins: acylation step.

The acylation step of the catalytic mechanism of beta-lactamases and penicillin-binding proteins (PBPs) has been studied with various approaches. The methods applied range from molecular dynamics (MD) simulations to multiple titration calculations using the Poisson-Boltzmann approach to quantum mechanical (QM) methods. The mechanism of class A beta-lactamases was investigated in the greatest detail. Most approaches support the critical role of Glu-166 and hydrolytic water in the acylation step of the enzymatic catalysis in class A beta-lactamases. The details of the catalytic mechanism have been revealed by the QM approach, which clearly pointed out the critical role of Glu-166 acting as a general base in the acylation step with preferred substrates. Lys-73 shuffles a proton abstracted by Glu-166 O(epsilon ) to the beta-lactam nitrogen through Ser-130 hydroxyl. This proton is transferred from O(gamma) of the catalytic Ser-70 through the bridging hydrolytic water to Glu-166 O(epsilon ). Then the hydrogen is simultaneously passed through S(N)2 inversion mechanism at Lys-73 N(zeta) to Ser-130 O(gamma), which loses its proton to the beta-lactam nitrogen. The protonation of beta-lactam nitrogen proceeds with an immediate ring opening and collapse of the first tetrahedral species into an acyl-enzyme intermediate. However, the studies that considered the effect of solvation lower the barrier for the pathway, which utilizes Lys-73 as a general base, thus creating a possibility of multiple mechanisms for the acylation step in the class A beta-lactamases. These findings help explain the exceptional efficiency of these enzymes. They emphasize an important role of Glu-166, Lys-73, and Ser-130 for enzymatic catalysis and shed light on details of the acylation step of class A beta-lactamase mechanism. The acylation step for class C beta-lactamases and six classes of PBPs were also considered with continuum solvent models and MD simulations.

Simulating proteins at constant pH: An approach combining molecular dynamics and Monte Carlo simulation.

For the structure and function of proteins, the pH of the solution is one of the determining parameters. Current molecular dynamics (MD) simulations account for the solution pH only in a limited way by keeping each titratable site in a chosen protonation state. We present an algorithm that generates trajectories at a Boltzmann distributed ensemble of protonation states by a combination of MD and Monte Carlo (MC) simulation. The algorithm is useful for pH-dependent structural studies and to investigate in detail the titration behavior of proteins. The method is tested on the acidic residues of the protein hen egg white lysozyme. It is shown that small structural changes may have a big effect on the pK(A) values of titratable residues.

Conformational preferences of substituted prolines in the collagen triple helix.

Researchers have recently questioned the role hydroxylated prolines play in stabilizing the collagen triple helix. To address these issues, we have developed new molecular mechanics parameters for the simulation of peptides containing 4(R)-fluoroproline (Flp), 4(R)-hydroxyproline (Hyp), and 4(R)-aminoproline (Amp). Simulations of peptides based on these parameters can be used to determine the components that stabilize hydroxyproline over proline in the triple helix. The dihedrals F-C-C-N, O-C-C-N, and N-C-C-N were built using a N-beta-ethyl amide model. One nanosecond simulations were performed on the trimers [(Pro-Pro-Gly)(10)](3), [(Pro-Hyp-Gly)(10)](3), [(Pro-Amp-Gly)(10)](3), [(Pro-Amp(1+)-Gly)(10)](3), and [(Pro-Flp-Gly)(10)](3) in explicit solvent. The results of our simulations suggest that pyrrolidine ring conformation is mediated by the strength of the gauche effect and classical electrostatic interactions.

Theoretical and experimental studies of biotin analogues that bind almost as tightly to streptavidin as biotin.

We have used a newly developed qualitative computational approach, PROFEC (Pictorial Representation of Free Energy Changes), to visualize the areas of the ligand biotin where modifications of its structure might lead to tighter binding to the protein streptavidin. The PROFEC analysis, which includes protein flexibility and ligand solvation/desolvation, led to the suggestion that the pro-9R hydrogen atom of biotin, which is in alpha-position to the CO(2)(-) group, might be changed to a larger group and lead to better binding with streptavidin and avidin. Free energy calculations supported this suggestion and predicted that the methyl analogue should bind approximately 3 kcal/mol more tightly to streptavidin, with this difference coming exclusively from the relative desolvation free energy of the ligand. The PROFEC analysis further suggested little or no improvement for changing the pro-9S hydrogen atom to a methyl group, and great reduction in changing the ureido N-H groups to N-CH(3). Stimulated by these results, we synthesized 9R-methylbiotin and 9S-methylbiotin, and their binding free energies and enthalpies were measured for interaction with streptavidin and avidin, respectively. In contrast to the calculated results, experiments found both 9-methylbiotin isomers to bind more weakly to streptavidin than biotin. The calculated preference for the binding of the 9R- over the 9S-stereoisomer was observed. In addition, 9-methylbiotin is considerably less soluble in water than biotin, as predicted by the calculation, and the 9R isomer is, to our knowledge, thus far the tightest binding analogue of biotin to streptavidin. Subsequently, X-ray structures of the complexes between streptavidin and both 9R- and 9S-methylbiotin were determined, and the structures were consistent with those used in the free energy calculations. Thus, the reason for the discrepancy between the calculated and experimental binding free energy does not lie in unusual binding modes for the 9-methylbiotins.

Computational alanine scanning of the 1:1 human growth hormone-receptor complex.

The MM-PBSA (Molecular Mechanics-Poisson-Boltzmann surface area) method was applied to the human Growth Hormone (hGH) complexed with its receptor to assess both the validity and the limitations of the computational alanine scanning approach. A 400-ps dynamical trajectory of the fully solvated complex was simulated at 300 K in a 101 A x 81 A x 107 A water box using periodic boundary conditions. Long-range electrostatic interactions were treated with the particle mesh Ewald (PME) summation method. Equally spaced snapshots along the trajectory were chosen to compute the binding free energy using a continuum solvation model to calculate the electrostatic desolvation free energy and a solvent-accessible surface area approach to treat the nonpolar solvation free energy. Computational alanine scanning was performed on the same set of snapshots by mutating the residues in the structural epitope of the hormone and the receptor to alanine and recomputing the deltaGbinding. To further investigate a particular structure, a 200-ps dynamical trajectory of an R43A hormone-receptor complex was simulated. By postprocessing a single trajectory of the wild-type complex, the average unsigned error of our calculated deltadeltaGbinding is approximately1 kcal/mol for the alanine mutations of hydrophobic residues and polar/charged residues without buried salt bridges. When residues involved in buried salt bridges are mutated to alanine, it is demonstrated that a separate trajectory of the alanine mutant complex can lead to reasonable agreement with experimental results. Our approach can be extended to rapid screening of a variety of possible modifications to binding sites.

Molecular dynamics and free energy analyses of cathepsin D-inhibitor interactions: insight into structure-based ligand design.

In this study, we compare the calculated and experimental binding free energies for a combinatorial library of inhibitors of cathepsin D (CatD), an aspartyl protease. Using a molecular dynamics (MD)-based, continuum solvent method (MM-PBSA), we are able to reproduce the experimental binding affinity for a set of seven inhibitors with an average error of ca. 1 kcal/mol and a correlation coefficient of 0.98. By comparing the dynamical conformations of the inhibitors complexed with CatD with the initial conformations generated by CombiBuild (University of California, San Francisco, CA, 1995), we have found that the docking conformation observed in an X-ray structure of one peptide inhibitor bound to CatD (Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 6796-6800) is in good agreement with our MD simulation. However, the DOCK scoring function, based on intermolecular van der Waals and electrostatics, using a distance-dependent dielectric constant (J. Comput. Chem. 1992, 13, 505-524), poorly reproduces the trend of experimental binding affinity for these inhibitors. Finally, the use of PROFEC (J. Comput.-Aided Mol. Des. 1998, 12, 215-227) analysis allowed us to identify two possible substitutions to improve the binding of one of the better inhibitors to CatD. This study offers hope that current methods of estimating the free energy of binding are accurate enough to be used in a multistep virtual screening protocol.

A classical and ab initio study of the interaction of the myosin triphosphate binding domain with ATP.

We used classical molecular mechanics (MM) simulations and quantum mechanical (QM) structural relaxations to examine the active site of myosin when bound to ATP. Two conformations of myosin have been determined by x-ray crystallography. In one, there is no direct interaction between switch 2 and the nucleotide (open state). In the other (closed state), the universally conserved switch 2 glycine forms a hydrogen bond with a gamma-phosphate oxygen. MM simulations indicate that the two states are thermodynamically stable and allow us to investigate the extent to which the P-loop, switch 1, and switch 2 are involved in hydrolysis. We find that the open structure has a higher affinity for ATP than the closed structure, and that ATP is distorted toward a transition state by interactions with the protein. We also examine how the structure of the binding site changes with either MgATP or CaATP as the nucleotide in myosin in the open conformer. Our analyses suggest that higher CaATPase rates occur because the leaving phosphate (P(i)) group is more weakly bound and dissociation occurs faster. Finally, we validate the use of a particular formulation of a QM methodology (Car-Parrinello) to further refine the structures of the active site.