The Influence of Starting Coordinates in Free Energy Simulations of Ligand Binding to Dihydrofolate Reductase

Abstract

Molecular dynamics (MD) simulation combined with free energy perturbation (FEP) methods have been used to study the key structural differences and relative free energies for the binding of 6-methyl-N5-deazapterin (N8 protonated) and the 8-substituted compound, 6,8-dimethyl-N5-deazapterin (N3 protonated), to dihydrofolate reductase (DHFR). The free energy changes have been calculated using a variety of initial X-ray coordinates derived from bacterial and vertebrate (including human) DHFRs, and both with and without the reduced cofactor nicotinamide adenine dinucleotide (NADPH) bound. Given a sufficiently long simulation time for the FEP calculations (ca. 200 ps), all structures obtained after mutating 6,8-methyl-N5-deazapterin to 6-methyl-N5-deazapterin exhibited hydrogen bond formation between a backbone carbonyl group of DHFR and H(N8) of 6-methyl-N5-deazapterin, analogous to that found in the X-ray crystal structure of N5-deazafolate(N8 protonated) bound to human DHFR. However, both simulation and experiment suggest this additional H-bonding does not greatly enhance thermodynamic stability, with experiment indicating at most a factor of 2 difference in the relative affinities of the two ligand cations for vertebrate DHFR. Moreover, a binding differential of 10 in favour of the protonated 8-substituted compound is found experimentally for bacterial DHFR. The MD/FEP calculations suggest that the relative cost of ligand desolvation may largely cancel the lowering of free energy obtained in the active site, resulting in predicted binding differences within the range indicated by the vertebrate and bacterial DHFR experiments. However, the theoretical free energy changes could not be obtained with the accuracy required for the rationalization of the observed species dependence. While sampling difficulties are known to be inherent in MD simulation methodologies, these studies with several initial coordinate sets have demonstrated the contribution of coordinate choice to this problem. The results indicate that for demanding protein-ligand binding problems such as this one, the accuracy of the method may be no better than ± 2 kcal/mol.     

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.     

Theoretical study of the ligand-CYP2B4 complexes: effect of structure on binding free energies and heme spin state

The molecular origins of temperature-dependent ligand-binding affinities and ligand-induced heme spin state conversion have been investigated using free energy analysis and DFT calculations for substrates and inhibitors of cytochrome P450 2B4 (CYP2B4), employing models of CYP2B4 based on CYP2C5(3LVdH)/CYP2C9 crystal structures, and the results compared with experiment. DFT calculations indicate that large heme-ligand interactions (ca. -15 kcal/mol) are required for inducing a high to low spin heme transition, which is correlated with large molecular electrostatic potentials (approximately -45 kcal/mol) at the ligand heteroatom. While type II ligands often contain oxygen and nitrogen heteroatoms that ligate heme iron, DFT results indicate that BP and MF heme complexes, with weak substrate-heme interactions (ca. -2 kcal/mol), and modest MEPS minima (>-35 kcal/mol) are high spin. In contrast, heme complexes of the CYP2B4 inhibitor, 4PI, the product of benzphetamine metabolism, DMBP, and water are low spin, have substantial heme-ligand interaction energies (<-15 kcal/mol) and deep MEPS minima (<-45 kcal/mol) near their heteroatoms. MMPBSA analysis of MD trajectories were made to estimate binding free energies of these ligands at the heme binding site of CYP2B4. In order to initially assess the realism of this approach, the binding free energy of 4PI inhibitor was computed and found to be a reasonable agreement with experiment: -7.7 kcal/mol [-7.2 kcal/mol (experiment)]. BP was determined to be a good substrate [-6.3 kcal/mol (with heme-ligand water), -7.3 kcal/mol (without ligand water)/-5.8 kcal/mol (experiment)], whereas the binding of MF was negligible, with only marginal binding binding free energy of -1.7 kcal/mol with 2-MF bound [-3.8 kcal/mol (experiment)], both with and without retained heme-ligand water. Analysis of the free energy components reveal that hydrophobic/nonpolar contributions account for approximately 90% of the total binding free energy of these substrates and are the source of their differential and temperature-dependent CYP2B4 binding. The results indicate the underlying origins of the experimentally observed differential binding affinities of BP and MF, and indicate the plausibility of the use of models derived from moderate sequence identity templates in conjunction with approximate free energy methods in the estimation of ligand-P450 binding affinities.     

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.     

Calculation of the Affinity of the λ Repressor-Operator Complex Based on Free Energy Component Analysis

A calculation of the binding free energy of the u repressor-operator complex is described based on free energy component analysis. The calculations are based on a thermodynamic cycle of seven steps decomposed into a total of 24 individual components. The values of these terms are estimated using a combination of empirical potential functions from AMBER, generalized Born - solvent accessibility calculations, elementary statistical mechanics and semiempirical physicochemical properties. Two alternative approaches are compared, one based on the crystal structure of the complex and the other based on the molecular dynamics simulation of the u repressor-operator complex. The calculated affinity is m 19.7 kcal/mol from the crystal structure calculation and m 17.9 kcal/mol from the MD method. The corresponding experimental affinity of the complex is about m 12.6 kcal/mol, indicating reasonable agreement between theory and experiment, considering the approximations involved in the computational methodology. The results are analyzed in terms of contributions from electrostatics, van der Waal interactions, the hydrophobic effect and solvent release. The capabilities and limitations of free energy component methodology are assessed and discussed on the basis of these results.     

Exhaustive mutagenesis in silico: multicoordinate free energy calculations on proteins and peptides

We have extended and applied a multicoordinate free energy method, chemical Monte Carlo/Molecular Dynamics (CMC/MD), to calculate the relative free energies of different amino acid side-chains. CMC/MD allows the calculation of the relative free energies for many chemical species from a single free energy calculation. We have previously shown its utility in host:guest chemistry (Pitera and Kollman, J Am Chem Soc 1998;120:7557-7567)1 and ligand design (Eriksson et al., J Med Chem 1999;42:868-881)2, and here demonstrate its utility in calculations of amino acid properties and protein stability. We first study the relative solvation free energies of N-methylated and acetylated alanine, valine, and serine amino acids. With careful inclusion of rotameric states, internal energies, and both the solution and vacuum states of the calculation, we calculate relative solvation free energies in good agreement with thermodynamic integration (TI) calculations. Interestingly, we find that a significant amount of the unfavorable solvation of valine seen in prior work (Sun et al., J Am Chem Soc 1992;114:6798-6801)3 is caused by restraining the backbone in an extended conformation. In contrast, the solvation free energy of serine is calculated to be less favorable than expected from experiment, due to the formation of a favorable intramolecular hydrogen bond in the vacuum state. These monomer calculations emphasize the need to accurately consider all significant conformations of flexible molecules in free energy calculations. This development of the CMC/MD method paves the way for computations of protein stability analogous to the biochemical technique of "exhaustive mutagenesis." We have carried out just such a calculation at position 133 of T4 lysozyme, where we use CMC/MD to calculate the relative stability of eight different side-chain mutants in a single free energy calculation. Our T4 calculations show good agreement with the prior free energy calculations of Veenstra et al. (Prot Eng 1997;10:789-807)4 and excellent agreement with the experiments of Mendel et al. (Science 1992;256:1798-1802).     

Structure and thermodynamics of RNA-protein binding: using molecular dynamics and free energy analyses to calculate the free energies of binding and conformational change

An adaptive binding mechanism, requiring large conformational rearrangements, occurs commonly with many RNA-protein associations. To explore this process of reorganization, we have investigated the conformational change upon spliceosomal U1A-RNA binding with molecular dynamics (MD) simulations and free energy analyses. We computed the energetic cost of conformational change in U1A-hairpin and U1A-internal loop binding using a hybrid of molecular mechanics and continuum solvent methods. Encouragingly, in all four free energy comparisons (two slightly different proteins, two different RNAs), the free macromolecule was more stable than the bound form by the physically reasonable value of approximately 10 kcal/mol. We calculated the absolute binding free energies for both complexes to be in the same range as that found experimentally.     

Ranking ligand binding affinities with avidin: a molecular dynamics-based interaction energy study.

The binding of 14 biotin analogues to avidin is examined to evaluate the viability of calculating binding free energy based on molecular dynamics (MD) trajectories. Two approaches were investigated in this work. The first one uses the linear interaction energy approximation, while the other approach utilizes the interaction free energy. The results obtained from these two methods were found to correlate well with the experimental binding free energy data for 10 out of 14 ligands. For the other four ligands, both methods overestimate their binding strength by more than 7 kcal/mol. Free energy calculations using the thermodynamic integration method are employed to understand this overestimation. The effect of protein flexibility on binding free energy calculation and the effect of charged or neutral ligands on the calculated results are discussed. MD simulations are shown to be able to provide insight into the interactions occurring in the active site and the origins of variations in binding free energy.     

Relative binding free energy calculations of inhibitors to two mutants (Glu46----Ala/Gln) of ribonuclease T1 using molecular dynamics/free energy perturbation approaches.

We present free energy perturbation calculations on the complexes of Glu46----Ala46 (E46A) and Glu46----Gln46 (E46Q) mutants of ribonuclease T1 (RNaseT1) with inhibitors 2'-guanosine monophosphate (GMP) and 2'-adenosine monophosphate (AMP) by a thermodynamic perturbation method implemented with molecular dynamics (MD). Using the available crystal structure of the RNaseT1-GMP complex, the structures of E46A-GMP and E46Q-GMP were model built and equilibrated with MD simulations. The structures of E46A-AMP and E46Q-AMP were obtained as a final structure of the GMP----AMP perturbation calculation respectively. The calculated difference in the free energy of binding (delta delta Gbind) was 0.31 kcal/mol for the E46A system and -1.04 kcal/mol for the E46Q system. The resultant free energies are much smaller than the experimental and calculated value of approximately 3 kcal/mol for the native RNaseT1, which suggests that both mutants have greater relative adenine affinities than native RNaseT1. Especially E46Q is calculated to have a larger affinity for adenine than guanine, as we suggested previously from the calculation on the native RNaseT1. Thus, the molecular dynamics/free energy perturbation method may be helpful in protein engineering, directed toward increasing or changing the substrate specificity of enzymes.     

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.     

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.     

Accurate PDZ/Peptide Binding Specificity with Additive and Polarizable Free Energy Simulations

PDZ domains contain 80–100 amino acids and bind short C-terminal sequences of target proteins. Their specificity is essential for cellular signaling pathways. We studied the binding of the Tiam1 PDZ domain to peptides derived from the C-termini of its Syndecan-1 and Caspr4 targets. We used free energy perturbation (FEP) to characterize the binding energetics of one wild-type and 17 mutant complexes by simulating 21 alchemical transformations between pairs of complexes. Thirteen complexes had known experimental affinities. FEP is a powerful tool to understand protein/ligand binding. It depends, however, on the accuracy of molecular dynamics force fields and conformational sampling. Both aspects require continued testing, especially for ionic mutations. For six mutations that did not modify the net charge, we obtained excellent agreement with experiment using the additive, AMBER ff99SB force field, with a root mean square deviation (RMSD) of 0.37 kcal/mol. For six ionic mutations that modified the net charge, agreement was also good, with one large error (3 kcal/mol) and an RMSD of 0.9 kcal/mol for the other five. The large error arose from the overstabilization of a protein/peptide salt bridge by the additive force field. Four of the ionic mutations were also simulated with the polarizable Drude force field, which represents the first test of this force field for protein/ligand binding free energy changes. The large error was eliminated and the RMS error for the four mutations was reduced from 1.8 to 1.2 kcal/mol. The overall accuracy of FEP indicates it can be used to understand PDZ/peptide binding. Importantly, our results show that for ionic mutations in buried regions, electronic polarization plays a significant role.     

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.     

A computational analysis of the binding affinities of FKBP12 inhibitors using the MM-PB/SA method

The FK506-binding proteins have been targets of pharmaceutical interests over years. We have studied the binding of a set of 12 nonimmunosuppressive small-molecule inhibitors to FKBP12 through molecular dynamics simulations. Each complex was subjected to 1-ns MD simulation conducted in an explicit solvent environment under constant temperature and pressure. The binding free energy of each complex was then computed by the MM-PB/SA method in the AMBER program. Our MM-PB/SA computation produced a good correlation between the experimentally determined and the computed binding free energies with a correlation coefficient (R(2)) of 0.93 and a standard deviation as low as 0.30 kcal/mol. The vibrational entropy term given by the normal mode analysis was found to be helpful for achieving this correlation. Moreover, an adjustment to one weight factor in the PB/SA model was essential to correct the absolute values of the final binding free energies to a reasonable range. A head-to-head comparison of our MM-PB/SA model with a previously reported Linear Response Approximation (LRA) model suggested that the MM-PB/SA method is more robust in binding affinity prediction for this class of compounds.     

Investigation of glucose binding sites on insulin

Possible insulin binding sites for D-glucose have been investigated theoretically by docking and molecular dynamics (MD) simulations. Two different docking programs for small molecules were used; Multiple Copy Simultaneous Search (MCSS) and Solvation Energy for Exhaustive Docking (SEED) programs. The configurations resulting from the MCSS search were evaluated with a scoring function developed to estimate the binding free energy. SEED calculations were performed using various values for the dielectric constant of the solute. It is found that scores emphasizing non-polar interactions gave a preferential binding site in agreement with that inferred from recent fluorescence and NMR NOESY experiments. The calculated binding affinity of -1.4 to -3.5 kcal/mol is within the measured range of -2.0 +/- 0.5 kcal/mol. The validity of the binding site is suggested by the dynamical stability of the bound glucose when examined with MD simulations with explicit solvent. Alternative binding sites were found in the simulations and their relative stabilities were estimated. The motions of the bound glucose during molecular dynamics simulations are correlated with the motions of the insulin side chains that are in contact with it and with larger scale insulin motions. These results raise the question of whether glucose binding to insulin could play a role in its activity. The results establish the complementarity of molecular dynamics simulations and normal mode analyses with the search for binding sites proposed with small molecule docking programs.     

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.     

Rapid protein-ligand docking using soft modes from molecular dynamics simulations to account for protein deformability: binding of FK506 to FKBP

Most current docking methods to identify possible ligands and putative binding sites on a receptor molecule assume a rigid receptor structure to allow virtual screening of large ligand databases. However, binding of a ligand can lead to changes in the receptor protein conformation that are sterically necessary to accommodate a bound ligand. An approach is presented that allows relaxation of the protein conformation in precalculated soft flexible degrees of freedom during ligand-receptor docking. For the immunosuppressant FK506-binding protein FKBP, the soft flexible modes are extracted as principal components of motion from a molecular dynamics simulation. A simple penalty function for deformations in the soft flexible mode is used to limit receptor protein deformations during docking that avoids a costly recalculation of the receptor energy by summing over all receptor atom pairs at each step. Rigid docking of the FK506 ligand binding to an unbound FKBP conformation failed to identify a geometry close to experiment as favorable binding site. In contrast, inclusion of the flexible soft modes during systematic docking runs selected a binding geometry close to experiment as lowest energy conformation. This has been achieved at a modest increase of computational cost compared to rigid docking. The approach could provide a computationally efficient way to approximately account for receptor flexibility during docking of large numbers of putative ligands and putative docking geometries.