Estimation of the binding affinities of FKBP12 inhibitors using a linear response method.

A series of non-immunosuppressive inhibitors of FK506 binding protein (FKBP12) are investigated using Monte Carlo statistical mechanics simulations. These small molecules may serve as scaffolds for chemical inducers of protein dimerization, and have recently been found to have FKBP12-dependent neurotrophic activity. A linear response model was developed for estimation of absolute binding free energies based on changes in electrostatic and van der Waals energies and solvent-accessible surface areas, which are accumulated during simulations of bound and unbound ligands. With average errors of 0.5 kcal/mol, this method provides a relatively rapid way to screen the binding of ligands while retaining the structural information content of more rigorous free energy calculations.     

Mapping heat exchange in an allosteric protein

Nicotinic acetylcholine receptors (AChRs) are synaptic ion channels that spontaneously isomerize (i.e., gate) between resting and active conformations. We used single-molecule electrophysiology to measure the temperature dependencies of mouse neuromuscular AChR gating rate and equilibrium constants. From these we estimated free energy, enthalpy, and entropy changes caused by mutations of amino acids located between the transmitter binding sites and the middle of the membrane domain. The range of equilibrium enthalpy change (13.4 kcal/mol) was larger than for free energy change (5.5 kcal/mol at 25°C). For two residues, the slope of the rate-equilibrium free energy relationship (Φ) was approximately constant with temperature. Mutant cycle analysis showed that both free energies and enthalpies are additive for energetically independent mutations. We hypothesize that changes in energy associated with changes in structure mainly occur close to the site of the mutation, and, hence, that it is possible to make a residue-by-residue map of heat exchange in the AChR gating isomerization. The structural correlates of enthalpy changes are discussed for 12 different mutations in the protein.     

Machine Learning Estimates of Natural Product Conformational Energies

Machine learning has been used for estimation of potential energy surfaces to speed up molecular dynamics simulations of small systems. We demonstrate that this approach is feasible for significantly larger, structurally complex molecules, taking the natural product Archazolid A, a potent inhibitor of vacuolar-type ATPase, from the myxobacterium Archangium gephyra as an example. Our model estimates energies of new conformations by exploiting information from previous calculations via Gaussian process regression. Predictive variance is used to assess whether a conformation is in the interpolation region, allowing a controlled trade-off between prediction accuracy and computational speed-up. For energies of relaxed conformations at the density functional level of theory (implicit solvent, DFT/BLYP-disp3/def2-TZVP), mean absolute errors of less than 1 kcal/mol were achieved. The study demonstrates that predictive machine learning models can be developed for structurally complex, pharmaceutically relevant compounds, potentially enabling considerable speed-ups in simulations of larger molecular structures.     

Interactions of Arg2 in the Mnt N-terminal arm with the central and flanking regions of the Mnt operator

Arg2, in the N-terminal arm of the Mnt repressor, plays an important role in determining operator-binding specificity. In the complex of the Mnt tetramer with the 21 base-pair mnt operator, there are four potential sites for Arg2 interactions, two in the central region of the operator and two on the outer flanks of the operator. Single-chain variants of the dimeric N-terminal domain of Mnt containing one Arg2 residue and one Lys2 or Met2 residue were constructed and interactions with operator DNA were probed using Fe. EDTA affinity cleavage. The results of these orientation studies show that the majority of the energetically significant interactions mediated by Arg2 occur in the central region of the mnt operator. The RK2, RA2, and RM2 mutations reduce the free energy of operator binding by 1.7 kcal/mol, 3.3 kcal/mol, and 4.9 kcal/mol, respectively. Double-mutant thermodynamic cycle analyses using the RA2, RM2, and operator variants also reveal interaction free energies between Arg2 and operator base-pairs 9, 10, 11, 12 and 13, which in aggregate account for most of the Arg2 contribution to operator binding.     

Predicting absolute ligand binding free energies to a simple model site

A central challenge in structure-based ligand design is the accurate prediction of binding free energies. Here we apply alchemical free energy calculations in explicit solvent to predict ligand binding in a model cavity in T4 lysozyme. Even in this simple site, there are challenges. We made systematic improvements, beginning with single poses from docking, then including multiple poses, additional protein conformational changes, and using an improved charge model. Computed absolute binding free energies had an RMS error of 1.9 kcal/mol relative to previously determined experimental values. In blind prospective tests, the methods correctly discriminated between several true ligands and decoys in a set of putative binders identified by docking. In these prospective tests, the RMS error in predicted binding free energies relative to those subsequently determined experimentally was only 0.6 kcal/mol. X-ray crystal structures of the new ligands bound in the cavity corresponded closely to predictions from the free energy calculations, but sometimes differed from those predicted by docking. Finally, we examined the impact of holding the protein rigid, as in docking, with a view to learning how approximations made in docking affect accuracy and how they may be improved.     

Molecular dynamics simulations of class C beta-lactamase from Citrobacter freundii: insights into the base catalyst for acylation

Herein, we present results from molecular dynamics (MD) simulations of the class C beta-lactamase from Citrobacter freundii and its Michaelis complex with aztreonam. Four different configurations of the active site were modeled in aqueous solution, and their relative stability was estimated by means of quantum mechanical energy calculations. For the free enzyme, the energetically most stable configurations present a neutral Lys67 residue or an anionic Tyr150 side chain. Our calculations predict that these two configurations are quite close in terms of free energy, the anionic Tyr150 state being favored by approximately 1 kcal/mol. In contrast, for the noncovalent complex formed between the C. freundii enzyme and aztreonam, the energetic analyses predict that the configuration with the neutral Lys67 residue is much more stable than the anionic Tyr150 one (approximately 20 kcal/mol). Moreover, the MD simulations reveal that the neutral Lys67 state results in a proper enzyme-aztreonam orientation for nucleophilic attack and in a very stable contact between the nucleophilic hydroxyl group of Ser64 and the neutral amino side chain of Lys67. Thus, both the computed free energies and the structural analyses support the assignation of Lys67 as the base catalyst for the acylation step in the native form of the C. freundii enzyme.     

Free energy of charges in solvated proteins: microscopic calculations using a reversible charging process

Evaluation of the free energy of ionization of acidic groups in proteins may be used as a powerful and general test case for determining the reliability of calculations of electrostatic energies in macromolecules. This work attacks this test case by using an adiabatic charging process that evaluates the changes in free energies associated with ionizing the acidic groups Asp-3 and Glu-7 in bovine pancreatic trypsin inhibitor and aspartic acid in solution. The results of these free energy calculations are very encouraging; the error range is about 1 kcal/mol for these free energy changes of about-70 kcal/mol. This indicates that we are finally approaching the stage of obtaining quantitative results in modeling the energetics of solvated proteins.     

Computational studies and peptidomimetic design for the human p53-MDM2 complex

The interaction between human p53 and MDM2 is a key event in controlling cell growth. Many studies have suggested that a p53 mimic would be sufficient to inhibit MDM2 to reduce cell growth in cancerous tissue. In order to design a potent p53 mimic, molecular dynamics (MD) simulations were used to examine the binding interface and the effect of mutating key residues in the human p53-MDM2 complex. The Generalized Born surface area (GBSA) method was used to estimate free energies of binding, and a computational alanine-scanning approach was used to calculate the relative effects in the free energy of binding for key mutations. Our calculations determine the free energy of binding for a model p53-MDM2 complex to be -7.4 kcal/mol, which is in very good agreement with the experimentally determined values (-6.6--8.8 kcal/mol). The alanine-scanning results are in good agreement with experimental data and calculations by other groups. We have used the information from our studies of human p53-MDM2 to design a beta-peptide mimic of p53. MD simulations of the mimic bound to MDM2 estimate a free energy of binding of -8.8 kcal/mol. We have also applied alanine scanning to the mimic-MDM2 complex and reveal which mutations are most likely to alter the binding affinity, possibly giving rise to escape mutants. The mimic was compared to nutlins, a new class of inhibitors that block the formation of the p53-MDM2 complex. There are interesting similarities between the nutlins and our mimic, and the differences point to ways that both inhibitors may be improved. Finally, an additional hydrophobic pocket is noted in the interior of MDM2. It may be possible to design new inhibitors to take advantage of that pocket.     

Free-energy simulations of the oxidation of c-terminal methionines in calmodulin

In the course of aging or under conditions of oxidative stress, methionine residues of calmodulin undergo oxidation, leading to loss of biological activity of the protein. We have performed free-energy simulations of the effects of C-terminal methionine side-chain oxidation on the properties of calmodulin. The simulation results indicate that oxidation should have a destabilizing effect on all three protein functional states: calcium free, calcium loaded, and with both calcium and target peptide bound. Because the different states are destabilized by different amounts, this leads to a more complex pattern in the observable effects on protein thermal stability, calcium affinity, and binding of a target peptide. The influence of oxidation on the free energy of CaM unfolding is estimated by comparing the free-energy cost of oxidizing a Met residue in a Gly-Met-Gly peptide and in the protein. The protein thermal stability of the oxidized forms is lowered by a moderate amount 1-3 kcal/mol, in qualitative agreement with experimental results of 0.3 kcal/mol. The calculated changes in affinity for calcium and for the target peptide show opposing trends. Oxidation at position 144 is predicted to enhance peptide binding and weaken calcium binding, whereas oxidation at 145 weakens peptide binding and enhances affinity for calcium. The lower affinity of Met 145-oxidized calmodulin toward the target peptide correlates with experimentally observed lowering of calmodulin-activated Ca-ATPase activity when oxidized calmodulin from aged rat brains is used. Thus, our simulations suggest that Met 145 is the oxidation site in the C-terminal fragment of calmodulin. The microscopic mechanism behind the calculated free energy changes appears to be a greater affinity for water of the oxidized Met side-chain relative to normal Met. Structures with Met exposed to solvent had consistently lower free energies than those with buried Met sidechains.     

Absolute binding free energy calculations using molecular dynamics simulations with restraining potentials

The absolute (standard) binding free energy of eight FK506-related ligands to FKBP12 is calculated using free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent. A number of features are implemented to improve the accuracy and enhance the convergence of the calculations. First, the absolute binding free energy is decomposed into sequential steps during which the ligand-surrounding interactions as well as various biasing potentials restraining the translation, orientation, and conformation of the ligand are turned "on" and "off." Second, sampling of the ligand conformation is enforced by a restraining potential based on the root mean-square deviation relative to the bound state conformation. The effect of all the restraining potentials is rigorously unbiased, and it is shown explicitly that the final results are independent of all artificial restraints. Third, the repulsive and dispersive free energy contribution arising from the Lennard-Jones interactions of the ligand with its surrounding (protein and solvent) is calculated using the Weeks-Chandler-Andersen separation. This separation also improves convergence of the FEP/MD calculations. Fourth, to decrease the computational cost, only a small number of atoms in the vicinity of the binding site are simulated explicitly, while all the influence of the remaining atoms is incorporated implicitly using the generalized solvent boundary potential (GSBP) method. With GSBP, the size of the simulated FKBP12/ligand systems is significantly reduced, from approximately 25,000 to 2500. The computations are very efficient and the statistical error is small ( approximately 1 kcal/mol). The calculated binding free energies are generally in good agreement with available experimental data and previous calculations (within approximately 2 kcal/mol). The present results indicate that a strategy based on FEP/MD simulations of a reduced GSBP atomic model sampled with conformational, translational, and orientational restraining potentials can be computationally inexpensive and accurate.     

Ab initio computational study of beta-cellobiose conformers using B3LYP/6-311++G**

The molecular structure of 27 conformers of beta-cellobiose were studied in vacuo through gradient geometry optimization using B3LYP density functionals and the 6-311++G** basis set. The conformationally dependent geometry changes and energies were explored as well as the hydrogen-bonding network. The lowest electronic energy structures found were not those suggested from available crystallographic and NMR solution data, where the glycosidic dihedral angles fall in the region (phi, psi) approximately (40 degrees, -20 degrees ). Rather, 'flipped' conformations in which the dihedral angles are in the range (phi, psi) approximately (180 degrees, 0 degrees ) are energetically more stable by approximately 2.5 kcal/mol over the 'experimentally accepted' structure. Further, when the vibrational free energy, deltaG, obtained from the calculated frequencies, is compared throughout the series, structures with (phi, psi) in the experimentally observed range still have higher free energy ( approximately 2.0 kcal/mol) than 'flipped' forms. The range of bridging dihedral angles of the 'normal' conformers, resulting from the variance in the phi dihedral is larger than that found in the 'flipped' forms. Due to this large flat energy surface for the normal conformations, we surmise that the summation of populations of these conformations will favor the 'normal' conformations, although evidence suggests that polar solvent effects may play the dominant role in providing stability for the 'normal' forms. Even though some empirical studies previously found the 'flipped' conformations to be lowest in energy, these studies have been generally discredited because they were in disagreement with experimental results. Most of the DFT/ab initio conformations reported here have not been reported previously in the ab initio literature, in part because the use of less rigorous theoretical methods, i.e. smaller basis sets, have given results in general agreement with experimental data, that is, they energetically favored the 'normal' forms. These are the first DFT/ab initio calculations at this level of theory, apparently because of the length and difficulty of carrying out optimizations at these high levels.     

Probing the effect of point mutations at protein-protein interfaces with free energy calculations

We have studied the effect of point mutations of the primary binding residue (P1) at the protein-protein interface in complexes of chymotrypsin and elastase with the third domain of the turkey ovomucoid inhibitor and in trypsin with the bovine pancreatic trypsin inhibitor, using molecular dynamics simulations combined with the linear interaction energy (LIE) approach. A total of 56 mutants have been constructed and docked into their host proteins. The free energy of binding could be reliably calculated for 52 of these mutants that could unambiguously be fitted into the binding sites. We find that the predicted binding free energies are in very good agreement with experimental data with mean unsigned errors between 0.50 and 1.03 kcal/mol. It is also evident that the standard LIE model used to study small drug-like ligand binding to proteins is not suitable for protein-protein interactions. Three different LIE models were therefore tested for each of the series of protein-protein complexes included, and the best models for each system turn out to be very similar. The difference in parameterization between small drug-like compounds and protein point mutations is attributed to the preorganization of the binding surface. Our results clearly demonstrate the potential of free energy calculations for probing the effect of point mutations at protein-protein interfaces and for exploring the principles of specificity of hot spots at the interface.     

Contributions to the binding free energy of ligands to avidin and streptavidin

The free energy of binding of a ligand to a macromolecule is here formally decomposed into the (effective) energy of interaction, reorganization energy of the ligand and the macromolecule, conformational entropy change of the ligand and the macromolecule, and translational and rotational entropy loss of the ligand. Molecular dynamics simulations with implicit solvation are used to evaluate these contributions in the binding of biotin, biotin analogs, and two peptides to avidin and streptavidin. We find that the largest contribution opposing binding is the protein reorganization energy, which is calculated to be from 10 to 30 kcal/mol for the ligands considered here. The ligand reorganization energy is also significant for flexible ligands. The translational/rotational entropy is 4.5-6 kcal/mol at 1 M standard state and room temperature. The calculated binding free energies are in the correct range, but the large statistical uncertainty in the protein reorganization energy precludes precise predictions. For some complexes, the simulations show multiple binding modes, different from the one observed in the crystal structure. This finding is probably due to deficiencies in the force field but may also reflect considerable ligand flexibility.     

B3LYP/6-311++G** study of alpha- and beta-D-glucopyranose and 1,5-anhydro-D-glucitol: 4C1 and 1C4 chairs, (3,O)B and B(3,O) boats, and skew-boat conformations

Geometry optimization, at the B3LYP/6-311++G** level of theory, was carried out on 4C1 and 1C4 chairs, (3,O)B and B(3,O) boats, and skew-boat conformations of alpha- and beta-D-glucopyranose. Similar calculations on 1,5-anhydro-D-glucitol allowed examination of the effect of removal of the 1-hydroxy group on the energy preference of the hydroxymethyl rotamers. Stable minimum energy boat conformers of glucose were found, as were stable skew boats, all having energies ranging from approximately 4-15 kcal/mol above the global energy 4C1 chair conformation. The 1C4 chair electronic energies were approximately 5-10 kcal/mol higher than the 4C1 chair, with the 1C4 alpha-anomers being lower in energy than the beta-anomers. Zero-point energy, enthalpy, entropy, and relative Gibbs free energies are reported at the harmonic level of theory. The alpha-anomer 4C1 chair conformations were found to be approximately 1 kcal/mol lower in electronic energy than the beta-anomers. The hydroxymethyl gt conformation was of lowest electronic energy for both the alpha- and beta-anomers. The glucose alpha/beta anomer ratio calculated from the relative free energies is 63/37%. From a numerical Hessian calculation, the tg conformations were found to be approximately 0.4-0.7 kcal/mol higher in relative free energy than the gg or gt conformers. Transition-state barriers to rotation about the C-5-C-6 bond were calculated for each glucose anomer with resulting barriers to rotation of approximately 3.7-5.8 kcal/mol. No energy barrier was found for the path between the alpha-gt and alpha-gg B(3,O) boat forms and the equivalent 4C1 chair conformations. The alpha-tg conformation has an energy minimum in the 1S3 twist form. Other boat and skew-boat forms are described. The beta-anomer boats retained their starting conformations, with the exception of the beta-tg-(3,O)B boat that moved to a skew form upon optimization.     

Calculation of cyclodextrin binding affinities: energy, entropy, and implications for drug design

The second generation Mining Minima method yields binding affinities accurate to within 0.8 kcal/mol for the associations of alpha-, beta-, and gamma-cyclodextrin with benzene, resorcinol, flurbiprofen, naproxen, and nabumetone. These calculations require hours to a day on a commodity computer. The calculations also indicate that the changes in configurational entropy upon binding oppose association by as much as 24 kcal/mol and result primarily from a narrowing of energy wells in the bound versus the free state, rather than from a drop in the number of distinct low-energy conformations on binding. Also, the configurational entropy is found to vary substantially among the bound conformations of a given cyclodextrin-guest complex. This result suggests that the configurational entropy must be accounted for to reliably rank docked conformations in both host-guest and ligand-protein complexes. In close analogy with the common experimental observation of entropy-enthalpy compensation, the computed entropy changes show a near-linear relationship with the changes in mean potential plus solvation energy.     

Analysis of the binding energies of testosterone, 5alpha-dihydrotestosterone,androstenedione and dehydroepiandrosterone sulfate with an antitestosterone antibody

Molecular dynamics simulations and molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) free energy calculations were used to study the binding of testosterone (TES), 5alpha-dihydrotestosterone (5ADHT), androstenedione (AND), and dehydroepiandrosterone sulfate (DHEAS) to the monoclonal antitestosterone antibody 3-C(4)F(5). The relative binding free energy of TES and AND was also calculated with free energy perturbation (FEP) simulations. The antibody 3-C(4)F(5) has a relatively high affinity (3 x 10(8) M(-1)) and on overall good binding profile for testosterone but its cross-reactivity with DHEAS has been the main reason for the failure to use this antibody in clinical immunoassays. The relative binding free energies obtained with the MM-PBSA method were 1.5 kcal/mol for 5ADHT, 3.8 kcal/mol for AND, and 4.3 kcal/mol for DHEAS, as compared to TES. When a water molecule of the ligand binding site, observed in the antibody-TES crystal structure, was explicitly included in MM-PBSA calculations, the relative binding energies were 3.4, 4.9, and 5.4 kcal/mol for 5ADHT, AND, and DHEAS, respectively. The calculated numbers are in correct order but larger than the corresponding experimental energies of 1.3, 1.5, and 2.6 kcal/mol, respectively. The fact that the MM-PBSA method reproduced the relative binding free energies of DHEAS, a steroid having a negatively charged sulfate group, and the neutrally charged TES, 5ADHT, and AND in satisfactory agreement with experiment shows the robustness of the method in predicting relative binding affinities. The 800-ps FEP simulations predicted that the antibody 3-C(4)F(5) binds TES 1.3 kcal/mol tighter than AND. Computational mutagenesis of selected amino acid residues of the ligand binding site revealed that the lower affinities of AND and DHEAS as compared to TES are due to a combined effect of several residues, each contributing a small fraction to the tighter binding of TES. An exception to this is Tyr99H, whose mutation to Ala lowered the binding of DHEAS 0.7 kcal/mol more than the binding of TES. This is probably due to the hydrogen bonding interaction formed between the OH group of Tyr99H and the sulfate group of DHEAS. Computational mutagensis data also showed that the affinity of the steroids to the antitestosterone antibody 3-C(4)F(5) would be enhanced if Trp47H were repositioned so that it would make more extensive contacts with the bound ligands. In addition, the binding of steroids to antitestosterone, antiprogesterone, and antiestradiol antibodies is discussed.     

Conformation of the Ras-binding domain of Raf studied by molecular dynamics and free energy simulations

Recognition of Ras by its downstream target Raf is mediated by a Ras-recognition region in the Ras-binding domain (RBD) of Raf. Residues 78–89 in this region occupy two different conformations in the ensemble of NMR solution structures of the RBD: a fully α-helical one, and one where 87–90 form a type IV β-turn. Molecular dynamics simulations of the RBD in solution were performed to explore the stability of these and other possible conformations of both the wild-type RBD and the R89K mutant, which does not bind Ras. The simulations sample a fully helical conformation for residues 78–89 similar to the NMR helical structures, a conformation where 85–89 form a 3₁₀-helical turn, and a conformation where 87–90 form a type I |iB-turn, whose free energies are all within 0.3 kcal/mol of each other. NOE patterns and H^α chemical shifts from the simulations are in reasonable agreement with experiment. The NMR turn structure is calculated to be 3 kcal/mol higher than the three above conformations. In a simulation with the same implicit solvent model used in the NMR structure generation, the turn conformation relaxes into the fully helical conformation, illustrating possible structural artifacts introduced by the implicit solvent model. With the Raf R89K mutant, simulations sample a fully helical and a turn conformation, the turn being 0.9 kcal/mol more stable. Thus, the mutation affects the population of RBD conformations, and this is expected to affect Ras binding. For example, if the fully helical conformation of residues 78–89 is required for binding, its free energy increase in R89K will increase the binding free energy by about 0.6 kcal/mol. Proteins 31:186–200, 1998. © 1998 Wiley-Liss, Inc.     

Probing carbohydrate product expulsion from a processive cellulase with multiple absolute binding free energy methods

Understanding the enzymatic mechanism that cellulases employ to degrade cellulose is critical to efforts to efficiently utilize plant biomass as a sustainable energy resource. A key component of cellulase action on cellulose is product inhibition from monosaccharide and disaccharides in the product site of cellulase tunnel. The absolute binding free energy of cellobiose and glucose to the product site of the catalytic tunnel of the Family 7 cellobiohydrolase (Cel7A) of Trichoderma reesei (Hypocrea jecorina) was calculated using two different approaches: steered molecular dynamics (SMD) simulations and alchemical free energy perturbation molecular dynamics (FEP/MD) simulations. For the SMD approach, three methods based on Jarzynski's equality were used to construct the potential of mean force from multiple pulling trajectories. The calculated binding free energies, -14.4 kcal/mol using SMD and -11.2 kcal/mol using FEP/MD, are in good qualitative agreement. Analysis of the SMD pulling trajectories suggests that several protein residues (Arg-251, Asp-259, Asp-262, Trp-376, and Tyr-381) play key roles in cellobiose and glucose binding to the catalytic tunnel. Five mutations (R251A, D259A, D262A, W376A, and Y381A) were made computationally to measure the changes in free energy during the product expulsion process. The absolute binding free energies of cellobiose to the catalytic tunnel of these five mutants are -13.1, -6.0, -11.5, -7.5, and -8.8 kcal/mol, respectively. The results demonstrated that all of the mutants tested can lower the binding free energy of cellobiose, which provides potential applications in engineering the enzyme to accelerate the product expulsion process and improve the efficiency of biomass conversion.