Free energy perturbation calculations on parallel computers: Demonstrations of scalable linear speedup

A coarse-grain parallel implementation of the free energy perturbation (FEP) module of the AMBER molecular dynamics program is described and then demonstrated using five different molecular systems. The difference in the free energy of (aqueous) solvation is calculated for two monovalent cations ΔΔG_aq(Li⁺ Δ Cs⁺), and for the zero-sum ethane-to-ethane′ perturbation ΔΔG_aq(CH₃? methyl? X → X? methyl? CH₃), where X is a ghost methyl. The difference in binding free energy for a docked HIV-1 protease inhibitor into its ethylene mimetic is examined by mutating its fifth peptide bond, ΔG(CO? NH → CH?CH). A potassium ion (K⁺) is driven outward from the center of mass of ionophore salinomycin (SAL^?) in a potential of mean force calculation ΔG_MeOH(SAL^? · K⁺) carried out in methanol solvent. Parallel speedup obtained is linearly proportional to the number of parallel processors applied. Finally, the difference in free energy of solvation of phenol versus benzene, ΔΔG_oct(phenol → benzene), is determined in water-saturated octanol and then expressed in terms of relative partition coefficients, Δ log(P_o/w). Because no interprocessor communication is required, this approach is scalable and applicable in general for any parallel architecture or network of machines. FEP calculations run on the nCUBE/2 using 50 or 100 parallel processors were completed in clock times equivalent to or twice as fast as a Cray Y-MP. The difficulty of ensuring adequate system equilibrium when agradual configurational reorientation follows the mutation of the Hamiltonian is discussed and analyzed. The results of a successful protocol for overcoming this equilibration problem are presented. The types of molecular perturbations for which this method is expected to perform most efficiently are described. © 1994 by John Wiley & Sons, Inc.     

Free energy perturbation and molecular dynamics calculations of copper binding to azurin

Free energy perturbation/molecular dynamics simulations have been carried out on copper/azurin systems calculating the binding affinities of copper (II) ion to azurin either in the native or in the unfolded state. In order to test the validity of the strategy adopted for the calculations and to establish what force field is suitable for these kinds of calculations, three different force fields, AMBER, CVFF, and CFF, have been alternatively used for the calculations and the results have been compared with experimental data obtained by spectroscopic titrations of copper (II)/azurin solutions and denaturation experiments. Our findings have pointed out that only CFF gives satisfactory results, thus providing a reliable tool for copper binding simulations in copper protein.     

Free energy perturbation approach to the critical assessment of selective cyclooxygenase-2 inhibitors

The discovery of selective cyclooxygenase-2 (COX-2) inhibitors represents a major achievement of the efforts over the past few decades to develop therapeutic treatments for inflammation. To gain insights into designing new COX-2-selective inhibitors, we address the energetic and structural basis for the selective inhibition of COX isozymes by means of a combined computational protocol involving docking experiment, force field design for the heme prothetic group, and free energy perturbation (FEP) simulation. We consider both COX-2- and COX-1-selective inhibitors taking the V523I mutant of COX-2 to be a relevant structural model for COX-1 as confirmed by a variety of experimental and theoretical evidences. For all COX-2-selective inhibitors under consideration, we find that free energies of binding become less favorable as the receptor changes from COX-2 to COX-1, due to the weakening and/or loss of hydrogen bond and hydrophobic interactions that stabilize the inhibitors in the COX-2 active site. On the other hand, COX-1-selective oxicam inhibitors gain extra stabilization energy with the change of residue 523 from valine to isoleucine because of the formations of new hydrogen bonds in the enzyme-inhibitor complexes. The utility of the combined computational approach, as a valuable tool for in silico screening of COX-2-selective inhibitors, is further exemplified by identifying the physicochemical origins of the enantiospecific selective inhibition of COX-2 by -substituted indomethacin ethanolamide inhibitors.     

Prediction of the potency of inhibitors of adenosine deaminase by QM/MM calculations

QM/MM calculations show that the potency of a range of inhibitors of adenosine deaminase correlates with the relative stability of the reaction intermediate at the active site, rather than with the inhibitor binding energy.     

Predicting relative binding affinities of non-peptide HIV protease inhibitors with free energy perturbation calculations

The relative binding free energies in HIV protease of haloperidol thioketal (THK) and three of its derivatives were examined with free energy calculations. THK is a weak inhibitor (IC50 = 15 M) for which two cocrystal structures with HIV type 1 proteases have been solved [Rutenber, E. et al., J. Biol. Chem., 268 (1993) 15343]. A THK derivative with a phenyl group on C2 of the piperidine ring was expected to be a poor inhibitor based on experiments with haloperidol ketal and its 2- phenyl derivative (Caldera, P., personal communication). Our calculations predict that a 5-phenyl THK derivative, suggested based on examination of the crystal structure, will bind significantly better than THK. Although there are large error bars as estimated from hysteresis, the calculations predict that the 5-phenyl substituent is clearly favored over the 2-phenyl derivative as well as the parent compound. The unfavorable free energies of solvation of both phenyl THK derivatives relative to the parent compound contributed to their predicted binding free energies. In a third simulation, the change in binding free energy for 5-benzyl THK relative to THK was calculated. Although this derivative has a lower free energy in the protein, its decreased free energy of solvation increases the predicted G(bind) to the same range as that of the 2-phenyl derivative.     

Free energy perturbation simulations of cation binding to valinomycin

Experimental values of the free energies of cation binding to the cyclic depsipeptide molecule, valinomycin, obtained from Pedersen-type salt extraction measurements, provide data against which it is possible to test the adequacy of the procedures and force fields of the molecular dynamics algorithms, MOLARIS and GROMOS. These data are then used to assess appropriate values for the partial charges of the ester carbonyl oxygen and carbon. Valinomycin was chosen because it has only one kind of ion-binding ligand and because the cation is sufficiently enfolded by the molecule in the ion-complexes that the overall size and shape of the complex is virtually the same regardless of the species of cation bound. For such an isosteric complex, the experimentally measured selectivities are sufficiently similar in a wide variety of solvent environments that thedifferences in free energies measured between the different ion-valinomycin complexes by two-phase salt extraction experiments into dichloromethane can be taken as equivalent to the differences in free energies in vacuo. Thesedifferences were therefore compared with those computed for ion-valinomycin complexationin vacuo by Free Energy Perturbation/Molecular Dynamics (FEP/MD) simulations using the MOLARIS and GROMOS programs. Starting with a set of Lennard-Jones 6–12 parameters for the monovalent cations assessed for aqueous solution we explored the effect of varying the partial charges of the ester carbonyl ligands on binding free energy differences (i.e. the selectivity) among Na, K, Rb, and Cs. The computed selectivity was found to depend strongly on the value of partial charge, following a typical Eisenman Selectivity Pattern in which the correct selectivity sequence and magnitude occurred only over a very narrow range of partial charge (around 0.33 and 0.6 for the standard carbonyls of MOLARIS and GROMOS, respectively). Using MOLARIS we explored the effect of varying the size of the ester carbonyl ligands by comparing the standard carbonyl of MOLARIS with the somewhat smaller carbonyls of GROMOS and found an equally satisfactory ability to reproduce the experimental data with a partial charge value of 0.41. These results validate the use of both the MOLARIS and GROMOS force fields as starting points for quantitative calculations of ion-binding in more complex molecules (e.g., ion-binding sites and channels in proteins).This paper is dedicated to the memory of the late Dr C. J. Pedersen.     

Free energy calculations on the two drug binding sites in the M2 proton channel

Two alternative binding sites of adamantane-type drugs in the influenza A M2 channel have been suggested, one with the drug binding inside the channel pore and the other with four drug molecule S-binding to the C-terminal surface of the transmembrane domain. Recent computational and experimental studies have suggested that the pore binding site is more energetically favorable but the external surface binding site may also exist. Nonetheless, which drug binding site leads to channel inhibition in vivo and how drug-resistant mutations affect these sites are not completely understood. We applied molecular dynamics simulations and potential of mean force calculations to examine the structures and the free energies associated with these putative drug binding sites in an M2-lipid bilayer system. We found that, at biological pH (~7.4), the pore binding site is more thermodynamically favorable than the surface binding site by ~7 kcal/mol and, hence, would lead to more stable drug binding and channel inhibition. This result is in excellent agreement with several recent studies. More importantly, a novel finding of ours is that binding to the channel pore requires overcoming a much higher energy barrier of ~10 kcal/mol than binding to the C-terminal channel surface, indicating that the latter site is more kinetically favorable. Our study is the first computational work that provides both kinetic and thermodynamic energy information on these drug binding sites. Our results provide a theoretical framework to interpret and reconcile existing and often conflicting results regarding these two binding sites, thus helping to expand our understanding of M2-drug binding, and may help guide the design and screening of novel drugs to combat the virus.     

Synthesis of 5'-methylthio coformycins: specific inhibitors for malarial adenosine deaminase

Transition state theory suggests that enzymatic rate acceleration (kcat/knon) is related to the stabilization of the transition state for a given reaction. Chemically stable analogues of a transition state complex are predicted to convert catalytic energy into binding energy. Because transition state stabilization is a function of catalytic efficiency, differences in substrate specificity can be exploited in the design of tight-binding transition state analogue inhibitors. Coformycin and 2'-deoxycoformycin are natural product transition state analogue inhibitors of adenosine deaminases (ADAs). These compounds mimic the tetrahedral geometry of the ADA transition state and bind with picomolar dissociation constants to enzymes from bovine, human, and protozoan sources. The purine salvage pathway in malaria parasites is unique in that Plasmodium falciparum ADA (PfADA) catalyzes the deamination of both adenosine and 5'-methylthioadenosine. In contrast, neither human adenosine deaminase (HsADA) nor the bovine enzyme (BtADA) can deaminate 5'-methylthioadenosine. 5'-Methylthiocoformycin and 5'-methylthio-2'-deoxycoformycin were synthesized to be specific transition state mimics of the P. falciparum enzyme. These analogues inhibited PfADA with dissociation constants of 430 and 790 pM, respectively. Remarkably, they gave no detectable inhibition of the human and bovine enzymes. Adenosine deamination is involved in the essential pathway of purine salvage in P. falciparum, and prior studies have shown that inhibition of purine salvage results in parasite death. Inhibitors of HsADA are known to be toxic to humans, and the availability of parasite-specific ADA inhibitors may prevent this side-effect. The potent and P. falciparum-specific inhibitors described here have potential for development as antimalarials without inhibition of host ADA.     

Alternative approaches to potential of mean force calculations: Free energy perturbation versus thermodynamic integration. Case study of some representative nonpolar interactions

We investigated the convergence behavior of potential of mean force (PMF) calculations using free energy perturbation (FEP), thermodynamic integration (TI), and “slow growth” (SG) techniques. The critical comparison of these alternative approaches is illustrated by the study of three different systems: two tagged argon atoms in a periodic box of argon, two methane molecules, and two benzene molecules maintained in a “T-shaped” conformation, both dimers embedded in a periodic box of water. The complete PMF simulations were carried out considering several protocols, in which the number of intermediate “λ” states, together with the amount of sampling per individual state, were varied. In most cases, as much as 1 ns of molecular dynamics (MD) sampling was used to derive each free energy profile. For the different systems examined, we find that FEP and TI unquestionably constitute robust computational methods leading to results of comparable accuracy. We also show that proper convergence of the free energy calculations, and further quantitative interpretation of the PMFs, requires total simulation times much higher than has been hitherto estimated. In some circumstances, the free energy profiles derived from FEP calculations tend to be slightly poorer than those obtained with TI, as a probable consequence of the greater sensitivity of FEP to the window spacing δλ. In the context of TI, and to a lesser extent FEP, simulations, it appears preferable to employ a limited number of “λ” points of the integrand involving extensive sampling, rather than numerous points with fewer samplings. Finally, we note that, at least in the case of nonpolar interactions, PMFs of reasonable quality can be generated using SG, and at a substantially lower cost than with either FEP or TI. © 1996 by John Wiley & Sons, Inc.     

Nucleoside free energy perturbation calculations: Mutation of purine-to-pyrimidine and pyrimidine-to-purine nucleosides

Free energy perturbation calculations were conducted on the mutations of pyrimidine-to-pyrimidine, purine-to-purine, purine-to-pyrimidine, and pyrimidine-to-purine nucleosides. The parameters and technique required for these perturbations is presented. Each of the four nucleosides in DNA were mutated into the other three nucleosides, and the calculated change in free energy for each of the 12 mutations is reported.     

Free energy perturbation study of water dimer dissociation kinetics

An efficient approach is described for using accurate ab initio calculations to determine the rates of elementary condensation and evaporation processes that lead to nucleation of aqueous aerosols. The feasibility of the method is demonstrated in an application to evaporation rates of water dimer at 230 K. The method, known as ABC-FEP (ab initio/classical free energy perturbation), begins with a calculation of the potential of mean force for the dissociation (evaporation) of small water clusters using a molecular dynamics (MD) simulation with a model potential. The free energy perturbation is used to calculate how changing from the model potential to a potential calculated from ab initio methods would alter the potential of mean force. The difference in free energy is the Boltzmann-weighted average of the difference between the ab initio and classical potential energies, with the average taken over a sample of configurations from the MD simulation. In principle, the method does not require a highly accurate model potential, though more accurate potentials require fewer configurations to achieve a small sampling error in the free energy perturbation step. To test the feasibility of obtaining accurate potentials of mean force from ab initio calculations at a modest number of configurations, the free energy perturbation method has been used to correct the errors when some standard models for bulk water (SPC, TIP4P, and TIP4PFQ) are applied to water dimer. To allow a thorough exploration of sampling issues, a highly accurate fit to results of accurate ab initio calculations, known as SAPT-5s, as been used a proxy for the ab initio calculations. It is shown that accurate values for a point on the potential of mean force can be obtained from any of the water models using ab initio calculations at only 50 configurations. Thus, this method allows accurate simulations of small clusters without the need to develop water models specifically for clusters.     

Free energy calculations on stiff chain constituents of polysaccharide gels

Summary The free energy, as a function of end-to-end length, is estimated for some of the polysaccharide constituents of pectate and alginate gels, for chains with degrees of polymerisation ranging from 10 to 1000. The energy is obtained from the distribution of end-to-end lengths estimated by a computer simulation method. Chains are generated in free space assuming the probability of a particular orientation occurring between adjacent monomers is related to the appropriate dimer energy by a Boltzmann function not; but, the dimer energy having been previously calculated by well established procedures.It is found that the 1–4 linked polymers of -D-galacturonic and -L-guluronic acids are appreciably stiffer than -D-mannuronic acid. However, for all three polymers, for degrees of polymerisation comparable to the shortest distances between crosslinks in alginate and pectate gels the chains are found to be constrained in an energy minimum close to the fully extended length. As the degree of polymerisation increases the shape of the free energy curve changes, though Gaussian behaviour is only observed in the case of polymannuronic chains with a degree of polymerisation of 1000.It is suggested that some of the properties of alginate and pectate gels can be understood in terms of the flexibilities of the polymer chains connecting adjacent crosslinks.     

Peptide mimetics as enzyme inhibitors: Use of free energy perturbation calculations to evaluate isosteric replacement for amide bonds in a potent HIV protease inhibitor

Summary We present the application of free energy perturbation theory/molecular dynamics to predict the consequence of replacing each of the seven peptide bonds in the potent HIV protease inhibitor JG365: ACE (acetyl)-Ser-Leu-Asn-HEA (hydroxyethylamine analog of Phe-Pro)-Ile-Val-NME (N-methyl) by ethylene or fluoroethylene isosteres. Replacing two of these bonds may well lead to significantly tighter binding; replacing two others is predicted to significantly diminish the binding affinity. Also, for three of the peptide bonds fluoroethylene replacements could lead to increased binding of free energies of the inhibitors. Our results should be considered as predictive since there are, as yet, no experimental results on such peptide replacements as enzyme inhibitors.     

Some remarks on the perturbation and variation–perturbation calculations of the free energy for quantum systems

A generalization of the Gibbs–Bogoliubov inequality F ? F₀ + 〈H ? H₀〉₀ for the free energy F is studied which leads to a variation principle for this quantity that may be of importance in certain computational applications to quantum systems. This approach is coupled with a study of the perturbation expansion of the free energy for a canonical ensemble with H = H₀ + λV in the general case when H₀ and V do not commute. The second- and high-order derivatives of the free energy with respect to the perturbation parameter λ are calculated. From the second-order term is finally obtained a second-order correction to the previous variational minimum for the free energy.     

Design and synthesis of new C-nucleosides as potential adenosine deaminase inhibitors

A number of new pyrazolo[3,4-c] and [4,3-b]pyridine C-nucleosides, which can be viewed as 4- or 6-deazaformycin analogues were synthesized and examined as potential adenosine deaminase (ADA) inhibitors. The compounds were prepared through the condensation of a suitably substituted, lithiated 2- or 4-methylpyridine with tri-O-benzyl-d-ribonolactone, followed by borohydride reduction of the resulting hemiacetals, intramolecular Mitsunobu cyclisation of the derived diols, formation of the pyrazolopyridine ring system and subsequent removal of the protecting groups. These derivatives were designed on the structural basis provided by docking simulations performed within the enzyme catalytic site, however they demonstrated weak ADA inhibitory activity. Theoretical calculations assisted in the interpretation of the obtained biological data, thus providing guidance for rational structural modifications within this molecular scaffold.     

Synthesis of fluorinated purine and 1-deazapurine glycosides as potential inhibitors of adenosine deaminase

The synthesis of 2- and 6-trifluoromethylated purines and 1-deazapurines was performed by formal [3 + 3]-cyclization reactions of 5-aminoimidazoles with a set of trifluoromethyl-substituted 1,3-CCC- and 1,3-CNC-dielectrophiles. The corresponding fluorinated nucleosides were synthesized by glycosylation of 9-unsubstituted purines and 1-deazapurines with peracetylated β-ribose, β-glucose, and rhamnose and subsequent deprotection. These scaffolds can be considered as potential inhibitors of adenosine deaminase (ADA) and inosine monophosphate dehydrogenase (IMPDH) enzymes.     

Application of SCF perturbation theory to molecular calculations

A method for solving Roothaan's molecular orbital equations by means of SCF perturbation theory is presented. An estimate of the accuracy of the third order expansion is made for the CNDO/2 approximation from a comparison of the results from direct calculations. It is found that the third order theory is sufficiently accurate for quantitative studies.