Are free energy calculations useful in practice? A comparison with rapid scoring functions for the p38 MAP kinase protein system.

Precise thermodynamic integration free energy simulations have been applied to a congeneric series of 16 inhibitors to the p38 MAP kinase protein for which the experimental binding data (IC(50)) is known. The relative free energy of binding for each compound has been determined. For comparison, the same series of compounds have also scored using the best rapid scoring functions used in database screening. From the results of these calculations, we find (1) that precise free energy simulations allow predictions that are reliable and in good agreement with experiment; (2) that predictions of lower reliability, but still in good qualitative agreement with experiment, can be obtained using the OWFEG free energy grid method, at a much lower computational cost; (3) and that other methods, not based on free energy simulations yield results in much poorer agreement with experiment. A new predictive index, which measures the reliability of a prediction method in the context of normal use, is defined and calculated for each scoring method. Predictive indices of 0.84, 0.56, 0.04, -0.05, and 0.25 are calculated for thermodynamic integration, OWFEG, ChemScore, PLPScore, and Dock Energy Score, respectively, where +1.0 is perfect correct prediction, -1.0 is perfect incorrect prediction, and 0.0 is random.

Improved scoring of ligand-protein interactions using OWFEG free energy grids.

A new approach to rapidly score protein-ligand interactions is tested on several protein-ligand systems. Results using this approach - the OWFEG free energy grid - are quite promising and are generally in better agreement with experiment (in some cases much better) than those obtained employing scoring techniques currently in wide use. The OWFEG free energy grid is generated from a one-window free energy perturbation MD simulation (Pearlman, D. A. J. Med. Chem. 1999, 42, 4313-4324). The OWFEG approach is applied to three protein systems: IMPDH, MAP kinase p38, and HIV-1 aspartyl protease. OWFEG scores are compared to experimental K(i) and IC50 data in each case. Empirical scoring functions applied to the same systems for comparison include ChemScore, Piecewise Linear Potential (PLP), and Dock energy score.

Free energy grids: a practical qualitative application of free energy perturbation to ligand design using the OWFEG method.

Traditional window-based free energy calculations can precisely determine the free energy corresponding to a molecular change of interest. However, calculations performed in this fashion are typically slow and resource-intensive, which renders them less than ideal for drug design. To circumvent this drawback, a new approximate free energy method, OWFEG, has been developed and tested. OWFEG replaces the exact free energy calculation for a single change with a set of approximate calculations for a grid of possible changes surrounding a molecule. One of the key features of OWFEG is that a floating independent reference frame (FIRF) is used, so that each grid point moves with the region of the molecule to which it is closest. In this way, this approach has been made applicable to flexible molecules. OWFEG is applied to two model systems and then to the FKBP-12.FK506 protein-ligand complex. On the basis of the results of these tests, this approximate method shows promise as a predictive tool for drug design.

Automated detection of problem restraints in NMR data sets using the FINGAR genetic algorithm method.

The recently described FINGAR genetic algorithm method for NMR refinement [D.A. Pearlman (1996) J. Biomol. NMR, 8, 67-76] has been extended so that it can be used to detect problem restraints in an NMR-derived set of data. A problem restraint is defined as a restraint in a generally well-behaved set where the associated target value is in error, due to inaccuracies in the data, misassignment, etc. The method described here, FINGAR.RWF, locates problem restraints by finding those restraints that, if removed from the data set, result in a disproportionate improvement in the scoring function. The method is applied to several test cases of simulated data, as well as to real data for the FK506 macrocycle, with excellent results.

Protein solution structure calculations in solution: solvated molecular dynamics refinement of calbindin D9k.

The three-dimensional solution structures of proteins determined with NMR-derived constraints are almost always calculated in vacuo. The solution structure of (Ca2+)2-calbindin D9k has been redetermined by new restrained molecular dynamics (MD) calculations that include Ca2+ ions and explicit solvent molecules. Four parallel sets of MD refinements were run to provide accurate comparisons of structures produced in vacuo, in vacuo with Ca2+ ions, and with two different protocols in a solvent bath with Ca2+ ions. The structural ensembles were analyzed in terms of structural definition, molecular energies, packing density, solvent-accessible surface, hydrogen bonds, and the coordination of calcium ions in the two binding loops. Refinement including Ca2+ ions and explicit solvent results in significant improvements in the precision and accuracy of the structure, particularly in the binding loops. These results are consistent with results previously obtained in free MD simulations of proteins in solution and show that the rMD refined NMR-derived solution structures of proteins, especially metalloproteins, can be significantly improved by these strategies.

CONCERTS: dynamic connection of fragments as an approach to de novo ligand design.

We have implemented and tested a new approach to de novo ligand design, CONCERTS (creation of novel compounds by evaluation of residues at target sites). In this method, each member of a user-defined set of fragments is allowed to move independently about a target active site during a molecular dynamics simulation. This allows the fragments to sample various low-energy orientations. When the geometry between proximal fragments is appropriate, bonds can be formed between the fragments. In this fashion, larger molecules can be built. The bonding arrangement can subsequently be changed-breaking bonds between chosen fragment pairs and forming them between other pairs-if the overall process creates lower energy molecules. We have tested this method with various mixes of fragments against the active sites of the FK506 binding protein (FKBP-12) and HIV-1 aspartyl protease. In several cases, CONCERTS suggests ligands which are in surprisingly good agreement with known inhibitors of these proteins.

FINGAR: A new genetic algorithm-based method for fitting NMR data.

A new NMR refinement method, FINGAR (FIt NMR using a Genetic AlgoRithm), has been developed, which allows one to determine a weighted set of structures that best fits measured NMR-derived data. This method shows appreciable advantages over commonly used refinement methods. FINGAR generates an ensemble of conformations whose average reproduces the experimental NMR-derived restraints. In addition, a statistical importance weight is assigned to each of the conformations in the ensemble. As a result, one is not limited to simply presenting an envelope of sampled conformers. Instead, one can subsequently focus on a select few conformers of high weight. This is critical, because many structural analyses depend on using discrete conformations, not simply averages or ensembles. The genetic algorithm used by FINGAR allows one to simultaneously and reliably fit against many restraints, and to generate solutions which include as many conformations with non-zero weights as are necessary to generate the best fit. An added benefit of FINGAR is that because the time-consuming step in this method needs only to be performed once, in the beginning of the first run, numerous FINGAR simulations can be performed rapidly.

Practical applications of time-averaged restrained molecular dynamics to ligand-receptor systems: FK506 bound to the Q50R,A95H,K98I triple mutant of FKBP-13.

The ability of time-averaged restrained molecular dynamics (TARMD) to escape local low-energy conformations and explore conformational space is compared with conventional simulated-annealing methods. Practical suggestions are offered for performing TARMD calculations with ligand-receptor systems, and are illustrated for the complex of the immunosuppressant FK506 bound to Q50R,A95H,K98I triple mutant FKBP-13. The structure of (13)C-labeled FK506 bound to triple-mutant FKBP-13 was determined using a set of 87 NOE distance restraints derived from HSQC-NOESY experiments. TARMD was found to be superior to conventional simulated-annealing methods, and produced structures that were conformationally similar to FK506 bound to wild-type FKBP-12. The individual and combined effects of varying the NOE restraint force constant, using an explicit model for the protein binding pocket, and starting the calculations from different ligand conformations were explored in detail.

Determination of the differential effects of hydrogen bonding and water release on the binding of FK506 to native and Tyr82-->Phe82 FKBP-12 proteins using free energy simulations.

We use the thermodynamic integration technique to calculate the free energy associated with the Tyr82-->Phe82 mutation (Y82F) in the protein FKBP-12, both free and bound to known inhibitor FK506 (tacrolimis). We find that the net difference in free energy for the two changes is 0.85 kcal/mol, with the binding of FK506 relatively more favorable for the native protein than the mutant. This net energy compares very favorably with the experimentally measured value of 0.60 kcal/mol. The results indicate that the relatively better binding of FK506 to the native protein is driven by the favorable entropy associated with the release of water molecules from the protein when the ligand binds. For a variety of reasons, modest size of the system, smallness of the change being examined, rapid convergence of the ensemble that needs to be determined and use of statistical estimates to control sampling, we have been able to carry out atypically reliable and reproducible free energy calculations for this protein system. Free energy changes for the two simulations (Y82F FKBP-12/FK506 and Y82F FKBP-12) have been calculated a total of eight times each, to compare a variety of different methodological choices and to ensure that the results are statistically significant. Detailed analysis of the free energy results has been carried out, and indicates that even when applicable, deconvolution of the total free energy into components can be very difficult, that the statistical error estimates can give a reasonable bound on the error in a simulation, and that one must be careful to use the same simulation protocol in all simulations being compared.

Solution structure of FK506 bound to the R42K, H87V double mutant of FKBP-12.

The binding of the FK506/FKBP-12 complex to calcineurin (CN), its putative target for immunosuppression, involves recognition of solvent-exposed regions of the ligand as well as FKBP-12 residues near the active site. The R42K, H87V double mutation of FKBP-12 decreases the CN affinity of the complex by 550-fold [Aldape, R. A., Futer, O., DeCenzo, M. T., Jarrett, B. P., Murcko, M. A., & Livingston, D. J. (1992) J. Biol. Chem. 267, 16029-16032]. This work reports the solution structure of 13C-labeled FK506 bound to R42K, H87V FKBP-12. Assignments and NOE measurements at three mixing times were made from inverse-detected 1H-13C NMR experiments. Structures were calculated by several different methods, including distance geometry, restrained molecular dynamics, and molecular dynamics with time-averaged restraints. The NMR structures of the ligand are very well defined by the NOE restraints and differ slightly from the X-ray structure in regions that are involved in crystal packing. Comparison with the NMR structure of FK506 bound to wild-type FKBP-12 reveals that the R42K, H87V mutation causes the ligand backbone near C16 to move by 2.5 to 4.5 A, reorients 15-MeO by 90 degrees, and shifts 13-MeO by approximately 1.5 A. FK506 appears to undergo a concerted, mutationally induced shift in the binding pocket, with the greatest changes occurring in the effector region of the drug. The altered effector conformation of mutant-bound FK506 may perturb interactions between the drug and CN, thus accounting for the effect of the double mutation upon the CN inhibitory activity of the complex.

How is an NMR structure best defined? An analysis of molecular dynamics distance-based approaches.

Model studies on the macrocyclic immunosuppressive agent FK506 challenge traditional approaches to defining a structure from data collected during a 2D NMR experiment. A variety of joint molecular dynamics/NMR-distance refinement methodologies are characterized. From the results it is clear that the traditional presentation of an NMR structure as a single representative minimized conformation or as a fairly tight envelope of conformers best meeting the imposed restraints can be misleading; a greater emphasis is required on dynamics and on the fact that an NMR structure represents a time average.

How well do time-averaged J-coupling restraints work?

A comparison is made of the consequences of using time-averaged and conventional vicinal (3)J-coupling restraints in molecular dynamics refinement of an adenosine nucleoside model system. The target values for the restraints are derived from a 3-ns unrestrained molecular dynamics simulation. A comparison of the results from the restrained refinements and the unrestrained trajectory reveals that while both restraint types (time-averaged and conventional) are capable of acceptably reproducing the averaged values of the restrained parameters, time-averaged J-coupling restraints allow a more realistic and thorough description of conformational fluctuations. The full description of conformational behavior for the sugar ring using time-averaged J-coupling restraints is in excellent agreement with the unrestrained results. J-coupling restraints can result in a localized 'heating effect' about the underlying torsion. This allows a restrained torsion to sample all low-energy rotomers separated by modest barriers in an appropriately weighted mixture that reproduces the J-restraint target value. This will generally be advantageous for experimentally derived data, though it can be misleading if all these low-energy rotomers did not contribute to the ensemble that yields the measured J-value. An analysis of how the force constant used in the restraint terms affects the refinement indicates that smaller force constants are to be preferred, and that constants in the range of K(j)≥0.4 kcal s(2)/mol are acceptably large to overcome the intrinsic preferences of the force field.

Free energy component analysis: application to the "Z-phobicity" of A-T base pairs.

We have carried out molecular dynamics/free energy perturbation calculations on the double helical hexamer d(CGCGCG)2 in both B and Z forms. The third C.G base pair was "mutated" to T-A in both B and Z-DNA. It is known experimentally that replacement of a C.G with a T-A base pair in an alternating CG sequence raises the energy of the Z form relative to the B form by approximately 1 kcal mole-1. We have carried out free energy component calculations to assess the reason for the "Z-phobicity" of T-A base pairs. There are two major contributions. The primary contribution is from the intra-base pair interactions of the mutated base pair itself which disfavor T-A relative to C.G in the Z form by congruent to 1.4 kcal mole-1. A secondary contribution of 0.4 kcal mole-1 arises because the two cytosines on the strand where G is mutated to A disfavor T-A relative to C.G in the Z form by 1.9 kcal mole-1, while the guanines flanking the C on the complementary strand stabilize the T-A base pair relative to the C.G in the Z form by 1.5 kcal mole-1. The effect of the phosphates, non-neighboring nucleotides and intramolecular energies of the base pair being mutated are all small and come close to canceling each other, leading to a net calculated free energy for T-A Z-phobicity of 1.7 kcal mole-1.

Comparison of solution structures of mutant bovine pancreatic trypsin inhibitor proteins using two-dimensional nuclear magnetic resonance.

Structural perturbations due to a series of mutations at the 30-51 disulfide bond of bovine pancreatic trypsin inhibitor have been explored using NMR. The mutants replaced cysteines at positions 30 and 51 by alanine at position 51 and alanine, threonine, or valine at position 30. Chemical shift changes occur in residues proximate to the site of mutation. NOE assignments were made using an automated procedure, NASIGN, which used information from the wild-type crystal structure. Intensity information was utilized by a distance geometry algorithm, VEMBED, to generate a series of structures for each protein. Statistical analyses of these structures indicated larger averaged structural perturbations than would be expected from crystallographic and other information. Constrained molecular dynamics refinement using AMBER at 900 K was useful in eliminating structural movements that were not a necessary consequence of the NMR data. In most cases, statistically significant movements are shown to be those greater than approximately 1 A. Such movements do not appear to occur between wild type and A30A51, a result confirmed by crystallography (Eigenbrot, C., Randal, M., & Kossiakoff, A.A., 1990, Protein Eng. 3, 591-598). Structural alterations in the T30A51 or V30A51 mutant proteins near the limits of detection occur in the beta-loop (residues 25-28) or C-terminal alpha-helix, respectively.

Solution structure of [d(GTATATAC)]2 via restrained molecular dynamics simulations with nuclear magnetic resonance constraints derived from relaxation matrix analysis of two-dimensional nuclear Overhauser effect experiments.

Two-dimensional nuclear Overhauser effect (2D NOE) spectra have been used as the experimental basis for determining the solution structure of the duplex [d(GTATATAC)]2 employing restrained molecular dynamics (rMD) simulations. The MARDIGRAS algorithm has been employed to construct a set of 233 interproton distance constraints via iterative complete relaxation matrix analysis utilizing the peak intensities from the 2D NOE spectra obtained for different mixing times and model structures. The upper and lower bounds for each of the constraints, defining size of a flat-well potential function term used in the rMD simulations, were conservatively chosen as the largest or smallest value calculated by MARDIGRAS. Three different starting models were utilized in several rMD calculations: energy-minimized A-DNA, B-DNA, and a structure containing wrinkled D-DNA in the interior. Considerable effort was made to define the appropriate force constants to be employed with the NOE terms in the AMBER force field, using as criteria the average constraints deviation, the constraints violation energy and the total energy. Of the 233 constraints, one was generated indirectly, but proved to be crucial in defining the structure: the cross-strand A5-H2 A5-H2 distance. As those two protons resonate isochronously for the self-complementary duplex, the distance cannot be determined directly. However, the general pattern of 2D NOE peak intensities, spin-lattice relaxation time (T1) values, and 31P nuclear magnetic resonance spectra lead to use of the A3-H2 A7-H2 distance for A5-H2 A5-H2 as well. Five rMD runs, with different random number seeds, were made for each of the three starting structures with the full distance constraint set. The average structure from all 15 runs and the five-structure averages from each starting structure were all quite similar. Two rMD runs for each starting structure were made with the A5-H2 A5-H2 constraint missing. The average of these six rMD runs revealed differences in structure, compared to that with the full set of constraints, primarily for the middle two base-pairs involving the missing cross-strand constraint but global deviations also were found. Conformational analysis of the resulting structures revealed that the inner four to six base-pairs differed in structure from the termini. Furthermore, an alternating structure was suggested with features alternating for the A-T and T-A steps.

Are time-averaged restraints necessary for nuclear magnetic resonance refinement? A model study for DNA.

A recently suggested method for refinement of structural data obtained from two-dimensional nuclear magnetic resonance experiments using molecular dynamics (MD) is explored. In this method, the time-averaged values of the appropriate internal co-ordinates of the molecule, calculated from the MD trajectory, are driven by restraints towards the experimental target values. This contrasts with most refinement procedures currently in use, where restraints are applied based on the instantaneous values of the appropriate co-ordinates. Both refinement methods are applied to the EcoRI restriction site DNA hexamer d(GAATTC)2, using target nuclear Overhauser enhancement distances derived from a one nanosecond unrestrained MD simulation of this structure. The resulting refined structures are compared to the results of the unrestrained MD trajectory, which serves as our "experimental" data. We show that although both methods can yield an average structure with the correct gross morphology, the new method allows both a much more realistic picture of inherent flexibility, and reproduces fine conformational detail better, such as sequence dependency. We also analyze the very long MD trajectory generated here (longer than any previously reported for a DNA oligomer), and find that significantly shorter simulations, typical of those frequently performed, may not yield acceptably reliable values for certain structural parameters.

The calculated free energy effects of 5-methyl cytosine on the B to Z transition in DNA.

We have examined the free energy effects of 5-methylation of cytosine on the B in equilibrium Z conformational equilibrium in DNA. Free energy differences were calculated using the free energy perturbation approach, which uses an easily derived equation from classical statistical mechanics to relate the free energy difference between two states to the ensemble average of the potential energy difference between the states. Calculations were carried both in explicit solvent and (for comparison) in vacuo. The free energy values obtained for the explicit solvent systems are total free energies, with contributions from all parts of the system (solvent + solute), and so are relevant to the B in equilibrium Z transitions observed under real (physiological) conditions. We calculate that in solution, methylation makes the B in equilibrium Z transition more favorable by about -0.4 kcal/mole base pair (bp) in free energy. This value compares well with approximate experimentally derived values of about -0.3 kcal/mole-bp. We also discuss a method for determining the free energy difference between conformational states poorly maintained by a potential energy model. Finally, the effects of methylation on the melting temperature of DNA are examined.

Why do A.T base pairs inhibit Z-DNA formation?

We have carried out free energy perturbation calculations on DNA double-stranded hexanucleotides. The sequence d(CGCGCG)2 has been "mutated" into d(CGTGCG).d(CGCACG) with the oligonucleotide in the A, B, and Z structural forms, both in vacuo and in aqueous solution. In addition, model free energy calculations have been carried out in which the electrostatic charges of the H-bonding groups of the bases in the major and minor grooves of the DNA are reduced to zero as a way of assessing the relative solvation effects of these groups in the different structural forms of DNA. Finally, energy component analyses have been carried out to assess the relative roles of different intranucleotide interactions on the B----Z equilibrium as a function of base sequence. In vacuo, the free energy for changing a G.C to an A.T base pair is largest in the Z conformation; in the A and B conformations, the free energy cost is approximately 2 kcal/mol lower (1 cal = 4.184 J). The results are similar when the simulations are run in explicit solvent: the change costs 3 kcal/mol more in the Z conformation than in the B form. These results are consistent with experimental data, where it is clear that A.T sequences are significantly more "Z-phobic" than G.C sequences. The calculations indicate that both intranucleotide and solvation interactions contribute to this Z-phobicity.

Atomic charges for DNA constituents derived from single-crystal X-ray diffraction data.

We have derived a complete set of atomic charges for DNA from very high resolution, low temperature, single-crystal X-ray diffraction data, collected for a variety of nucleosides and nucleotides: cytidine; deoxycytidine 5'-monophosphate; deoxythymidine; guanosine 5'-monophosphate; deoxyadenosine; adenosine. This set of charges represents the first experimentally based parameterization of an important term in the energy function used in most modeling of DNA. The resulting charges are in good agreement with chemical intuition and experimental observations. They also agree qualitatively with the theoretically derived values now commonly used, but numerous and significant quantitative differences are observed. Possible reasons for the quantitative disagreement are discussed. An averaged set of charges (derived from the experimental results), which can be used in DNA modeling calculations, is presented.

Effects of nucleotide bromination on the stabilities of Z-RNA and Z-DNA: a molecular mechanics/thermodynamic perturbation study.

The structures of ZI- and ZII-form RNA and DNA oligonucleotides were energy minimized in vacuum using the AMBER molecular mechanics force field. Alternating C-G sequences were studied containing either unmodified nucleotides, 8-bromoguanosine in place of all guanosine residues, 5-bromocytidine in place of all cytidine residues, or all modified residues. Some molecules were also energy minimized in the presence of H2O and cations. Free energy perturbation calculations were done in which G8 and C5 hydrogen atoms in one or two residues of Z-form RNAs and DNAs were replaced in a stepwise manner by bromines. Bromination had little effect on the structures of the energy-minimized molecules. Both the minimized molecular energies and the results of the perturbation calculations indicate that bromination of guanosine at C8 will stabilize the Z forms of RNA and DNA relative to the nonbrominated Z form, while bromination of cytidine at C5 stabilizes Z-DNA and destabilizes Z-RNA. These results are in agreement with experimental data. The destabilizing effect of br5C in Z-RNAs is apparently due to an unfavorable interaction between the negatively charged C5 bromine atom and the guanosine hydroxyl group. The vacuum-minimized energies of the ZII-form oligonucleotides are lower than those of the corresponding ZI-form molecules for both RNA and DNA. Previous x-ray diffraction, nmr, and molecular mechanics studies indicate that hydration effects may favor the ZI conformation over the ZII form in DNA. Molecular mechanics calculations show that the ZII-ZI energy differences for the RNAs are greater than three times those obtained for the DNAs. This is due to structurally reinforcing hydrogen-bonding interactions involving the hydroxyl groups in the ZII form, especially between the guanosine hydroxyl hydrogen atom and the 3'-adjacent phosphate oxygen. In addition, the cytidine hydroxyl oxygen forms a hydrogen bond with the 5'-adjacent guanosine amino group in the ZII-form molecule. Both of these interactions are less likely in the ZI-form molecule: the former due to the orientation of the GpC phosphate away from the guanosine ribose in the ZI form, and the latter apparently due to competitive hydrogen bonding of the cytidine 2'-hydroxyl hydrogen with the cytosine carbonyl oxygen in the ZI form. The hydrogen-bonding interaction between the cytidine hydroxyl oxygen and the 5'-adjacent guanosine amino group in Z-RNA twists the amino group out of the plane of the base. This may be responsible for differences in the CD and Raman spectra of Z-RNA and Z-DNA.