Assessment and acceleration of binding energy calculations for protein–ligand complexes by the fragment molecular orbital method

In the field of drug discovery, it is important to accurately predict the binding affinities between target proteins and drug applicant molecules. Many of the computational methods available for evaluating binding affinities have adopted molecular mechanics‐based force fields, although they cannot fully describe protein–ligand interactions. A noteworthy computational method in development involves large‐scale electronic structure calculations. Fragment molecular orbital (FMO) method, which is one of such large‐scale calculation techniques, is applied in this study for calculating the binding energies between proteins and ligands. By testing the effects of specific FMO calculation conditions (including fragmentation size, basis sets, electron correlation, exchange‐correlation functionals, and solvation effects) on the binding energies of the FK506‐binding protein and 10 ligand complex molecule, we have found that the standard FMO calculation condition, FMO2‐MP2/6‐31G(d), is suitable for evaluating the protein–ligand interactions. The correlation coefficient between the binding energies calculated with this FMO calculation condition and experimental values is determined to be R = 0.77. Based on these results, we also propose a practical scheme for predicting binding affinities by combining the FMO method with the quantitative structure–activity relationship (QSAR) model. The results of this combined method can be directly compared with experimental binding affinities. The FMO and QSAR combined scheme shows a higher correlation with experimental data (R = 0.91). Furthermore, we propose an acceleration scheme for the binding energy calculations using a multilayer FMO method focusing on the protein–ligand interaction distance. Our acceleration scheme, which uses FMO2‐HF/STO‐3G:MP2/6‐31G(d) at R_int = 7.0 Å, reduces computational costs, while maintaining accuracy in the evaluation of binding energy. © 2015 Wiley Periodicals, Inc.     

An analysis of the interactions between the Sem-5 SH3 domain and its ligands using molecular dynamics, free energy calculations, and sequence analysis.

The Src-homology-3 (SH3) domain of the Caenorhabditis elegans protein Sem-5 binds proline-rich sequences. It is reported that the SH3 domains broadly accept amide N-substituted residues instead of only recognizing prolines on the basis of side chain shape or rigidity. We have studied the interactions between Sem-5 and its ligands using molecular dynamics (MD), free energy calculations, and sequence analysis. Relative binding free energies, estimated by a method called MM/PBSA, between different substitutions at sites -1, 0, and +2 of the peptide are consistent with the experimental data. A new method to calculate atomic partial charges, AM1-BCC method, is also used in the binding free energy calculations for different N-substitutions at site -1. The results are very similar to those obtained from widely used RESP charges in the AMBER force field. AM1-BCC charges can be calculated more rapidly for any organic molecule than can the RESP charges. Therefore, their use can enable a broader and more efficient application of the MM/PBSA method in drug design. Examination of each component of the free energy leads to the construction of van der Waals interaction energy profiles for each ligand as well as for wild-type and mutant Sem-5 proteins. The profiles and free energy calculations indicate that the van der Waals interactions between the ligands and the receptor determine whether an N- or a Calpha-substituted residue is favored at each site. A VC value (defined as a product of the conservation percentage of each residue and its van der Waals interaction energy with the ligand) is used to identify several residues on the receptor that are critical for specificity and binding affinity. This VC value may have a potential use in identifying crucial residues for any ligand-protein or protein-protein system. Mutations at two of those crucial residues, N190 and N206, are examined. One mutation, N190I, is predicted to reduce the selectivity of the N-substituted residue at site -1 of the ligand and is shown to bind similarly with N- and Calpha-substituted residues at that site.     

Communication: Quantum polarized fluctuating charge model: a practical method to include ligand polarizability in biomolecular simulations

We present a simple and practical method to include ligand electronic polarization in molecular dynamics (MD) simulation of biomolecular systems. The method involves periodically spawning quantum mechanical (QM) electrostatic potential (ESP) calculations on an extra set of computer processors using molecular coordinate snapshots from a running parallel MD simulation. The QM ESPs are evaluated for the small-molecule ligand in the presence of the electric field induced by the protein, solvent, and ion charges within the MD snapshot. Partial charges on ligand atom centers are fit through the multi-conformer restrained electrostatic potential (RESP) fit method on several successive ESPs. The RESP method was selected since it produces charges consistent with the AMBER/GAFF force-field used in the simulations. The updated charges are introduced back into the running simulation when the next snapshot is saved. The result is a simulation whose ligand partial charges continuously respond in real-time to the short-term mean electrostatic field of the evolving environment without incurring additional wall-clock time. We show that (1) by incorporating the cost of polarization back into the potential energy of the MD simulation, the algorithm conserves energy when run in the microcanonical ensemble and (2) the mean solvation free energies for 15 neutral amino acid side chains calculated with the quantum polarized fluctuating charge method and thermodynamic integration agree better with experiment relative to the Amber fixed charge force-field.     

A numerically stable restrained electrostatic potential charge fitting method

Inspired by the idea of charge decomposition in calculation of the dipole preserving and polarization consistent charges (Zhang et al., J. Comput. Chem. 2011, 32, 2127), we have proposed a numerically stable restrained electrostatic potential (ESP)‐based charge fitting method for protein. The atomic charge is composed of two parts. The dominant part is fixed to a predefined value (e.g., AMBER charge), and the residual part is to be determined by restrained fitting to residual ESP on grid points around the molecule. Nonuniform weighting factors as a function of the dominant charge are assigned to the atoms. Because the residual part is several folds to several orders smaller than the dominant part, the impact of ill‐conditioning is alleviated. This charge fitting method can be used in quantum mechanical/molecular mechanical (QM/MM) simulations and similar studies, where QM calculated electronic properties are frequently mapped to partial atomic charges. © 2012 Wiley Periodicals, Inc.     

Force field development for cofactors in the photosystem II

We present a set of force field (FF) parameters compatible with the AMBER03 FF to describe five cofactors in photosystem II (PSII) of oxygenic photosynthetic organisms: plastoquinone‐9 (three redox forms), chlorophyll‐a, pheophytin‐a, heme‐b, and β‐carotene. The development of a reliable FF for these cofactors is an essential step for performing molecular dynamics simulations of PSII. Such simulations are important for the calculation of absorption spectrum and the further investigation of the electron and energy transfer processes. We have derived parameters for partial charges, bonds, angles, and dihedral‐angles from solid theoretical models using systematic quantum mechanics (QM) calculations. We have shown that the developed FF parameters are in good agreement with both ab initio QM and experimental structural data in small molecule crystals as well as protein complexes. © 2012 Wiley Periodicals, Inc.     

Extension of QM/MM docking and its applications to metalloproteins

To overcome the limitation of conventional docking methods which assume fixed charge model from force field parameters, combined quantum mechanics/molecular mechanics (QM/MM) method has been applied to docking as a variable charge model and shown to exhibit improvement on the docking accuracy over fixed charge based methods. However, it has also been shown that there are a number of examples for which adoption of variable‐charge model fails to reproduce the native binding modes. In particular, for metalloproteins, previously implemented method of QM/MM docking failed most often. This class of proteins has highly polarized binding sites at which high‐coordinate‐numbered metal ions reside. We extend the QM/MM docking method so that protein atoms surrounding the binding site along with metal ions are included as quantum region, as opposed to only ligand atoms. This extension facilitates the required scaling of partial charges on metal ions leading to prediction of correct binding modes in metalloproteins. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

The SAMPL4 hydration challenge: evaluation of partial charge sets with explicit-water molecular dynamics simulations

We used blind predictions of the 47 hydration free energies in the SAMPL4 challenge to test multiple partial charge models in the context of explicit solvent free energy simulations with the general AMBER force field. One of the partial charge models, IPolQ-Mod, is a fast continuum solvent-based implementation of the IPolQ approach. The AM1-BCC, restrained electrostatic potential (RESP) and IpolQ-Mod approaches all perform reasonably well (R² > 0.8), while VCharge, though faster, gives less accurate results (R² of 0.5). The AM1-BCC results are more accurate than those of RESP for tertiary amines and nitrates, but the overall difference in accuracy between these methods is not statistically significant. Interestingly, the IPolQ-Mod method is found to yield partial charges in very close agreement with RESP. This observation suggests that the success of RESP may be attributed to its fortuitously approximating the arguably more rigorous IPolQ approach.     

A three-point method for evaluations of AMBER force field parameters: an application to copper-based artificial nucleases

We present the theoretical evaluation of new AMBER force field parameters for 12 copper-based nucleases with bis(2-pyridylmethyl) amine, 2,2′-dipyridylamine, imidazole, N,N-bis(2-benzimidazolylmethyl) amine and their derivative ligands based on first-principles electronic structure calculations at the B3LYP level of theory. A three-point approach was developed to accurately and efficiently evaluate the force field parameters for the copper-based nucleases with the ligands. The protocol of RESP atomic charges has been used to calculate the atomic charge distributions of the studied copper-based nucleases. The evaluated force field parameters and RESP atomic charges have been successfully applied in the testing molecular mechanics calculations and molecular dynamics simulations on the nucleases and the nuclease–DNA complexes, respectively. It has been demonstrated that the developed force field parameters and atomic charges can consistently reproduce molecular geometries and conformations in the available X-ray crystal structures and can reasonably predict the interaction properties of the nucleases with DNA. The developed force field parameters in this work provide an extension of the AMBER force field for its application to computational modeling and simulations of the copper-based artificial nucleases associated with DNA.     

Binding structures of tri‐N‐acetyl‐β‐glucosamine in hen egg white lysozyme using molecular dynamics with a polarizable force field

Lysozyme is a well‐studied enzyme that hydrolyzes the β‐(1,4)‐glycosidic linkage of N‐acetyl‐β‐glucosamine (NAG)_n oligomers. The active site of hen egg‐white lysozyme (HEWL) is believed to consist of six subsites, A‐F that can accommodate six sugar residues. We present studies exploring the use of polarizable force fields in conjunction with all‐atom molecular dynamics (MD) simulations to analyze binding structures of complexes of lysozyme and NAG trisaccharide, (NAG)₃. MD trajectories are applied to analyze structures and conformation of the complex as well as protein–ligand interactions, including the hydrogen‐bonding network in the binding pocket. Two binding modes (ABC and BCD) of (NAG)₃ are investigated independently based on a fixed‐charge model and a polarizable model. We also apply molecular mechanics with generalized born and surface area (MM‐GBSA) methods based on MD using both nonpolarizable and polarizable force fields to compute binding free energies. We also study the correlation between root‐mean‐squared deviation and binding free energies of the wildtype and W62Y mutant; we find that for this prototypical system, approaches using the MD trajectories coupled with implicit solvent models are equivalent for polarizable and fixed‐charge models. © 2012 Wiley Periodicals, Inc.     

Dependency of ligand free energy landscapes on charge parameters and solvent models

We performed replica-exchange molecular dynamics (REMD) simulations of six ligands to examine the dependency of their free energy landscapes on charge parameters and solvent models. Six different charge parameter sets for each ligand were first generated by RESP and AM1-BCC methods using three different conformations independently. RESP charges showed some conformational dependency. On the other hand, AM1-BCC charges did not show conformational dependency and well reproduced the overall trend of RESP charges. The free energy landscapes obtained from the REMD simulations of ligands in vacuum, Generalized-Born (GB), and TIP3P solutions were then analyzed. We found that even small charge differences can produce qualitatively different landscapes in vacuum condition, but the differences tend to be much smaller under GB and TIP3P conditions. The simulations in the GB model well reproduced the landscapes in the TIP3P model using only a fraction of the computational cost. The protein-bound ligand conformations were rarely the global minimum states, but similar conformations were found to exist in aqueous solution without proteins in regions close to the global minimum, local minimum or intermediate states.     

Ab initio quantum mechanical study of the binding energies of human estrogen receptor alpha with its ligands: an application of fragment molecular orbital method

We have theoretically examined the relative binding affinities (RBA) of typical ligands, 17beta-estradiol (EST), 17alpha-estradiol (ESTA), genistein (GEN), raloxifene (RAL), 4-hydroxytamoxifen (OHT), tamoxifen (TAM), clomifene (CLO), 4-hydroxyclomifene (OHC), diethylstilbestrol (DES), bisphenol A (BISA), and bisphenol F (BISF), to the alpha-subtype of the human estrogen receptor ligand-binding domain (hERalpha LBD), by calculating their binding energies. The ab initio fragment molecular orbital (FMO) method, which we have recently proposed for the calculations of macromolecules such as proteins, was applied at the HF/STO-3G level. The receptor protein was primarily modeled by 50 amino acid residues surrounding the ligand. The number of atoms in these model complexes is about 850, including hydrogen atoms. For the complexes with EST, RAL, OHT, and DES, the binding energies were calculated again with the entire ERalphaLBD consisting of 241 residues or about 4000 atoms. No significant difference was found in the calculated binding energies between the model and the real protein complexes. This indicates that the binding between the protein and its ligands is well characterized by the model protein with the 50 residues. The calculated binding energies relative to EST were very well correlated with the experimental RBA (the correlation coefficient r=0.837) for the ligands studied in this work. We also found that the charge transfer between ER and ligands is significant on ER-ligand binding. To our knowledge, this is the first achievement of ab initio quantum mechanical calculations of large molecules such as the entire ERalphaLBD protein.     

'Inductive' charges on atoms in proteins: comparative docking with the extended steroid benchmark set and discovery of a novel SHBG ligand

We have developed a novel iterative approach for calculation of partial charges in proteins within the framework of the 'molecular capacitance' model. The method operates by an effective 'inductive' electronegativity scale derived from a number of the conventional charge systems including CHARMM, AMBER, MMFF, OPLS, and PEOE among others. Our novel 'inductive' electronegativity equalization procedure allows rapid and conformation sensitive computation of adequate partial charges in proteins. Accuracy of the 'inductive' values was confirmed by their correlation with DFT-computed partial charges in common amino acids. A comparative docking study with an extended steroid data set not only illustrated the adequacy of 'inductive' protein charges but also demonstrated their superior performance compared to several conventional protein charging systems. Subsequent docking with 'inductive' charges resulted in identification of five potential leads as human Sex Hormone Binding Globulin (SHBG) ligands from a commercial library of natural compounds. When the selected substances were evaluated for their ability to bind SHBG in vitro, three of them displaced testosterone from the SHBG steroid-binding site, and with one compound this was achieved at micromolar concentrations.     

Comparison of radii sets,entropy, QM methods,and sampling on MM‐PBSA,MM‐GBSA,and QM/MM‐GBSA ligand binding energies of F. tularensis enoyl‐ACP reductase (FabI)

To validate a method for predicting the binding affinities of FabI inhibitors, three implicit solvent methods, MM‐PBSA, MM‐GBSA, and QM/MM‐GBSA were carefully compared using 16 benzimidazole inhibitors in complex with Francisella tularensis FabI. The data suggests that the prediction results are sensitive to radii sets, GB methods, QM Hamiltonians, sampling protocols, and simulation length, if only one simulation trajectory is used for each ligand. In this case, QM/MM‐GBSA using 6 ns MD simulation trajectories together with GB^neck2, PM3, and the mbondi2 radii set, generate the closest agreement with experimental values (r² = 0.88). However, if the three implicit solvent methods are averaged from six 1 ns MD simulations for each ligand (called “multiple independent sampling”), the prediction results are relatively insensitive to all the tested parameters. Moreover, MM/GBSA together with GB^HCT and mbondi, using 600 frames extracted evenly from six 0.25 ns MD simulations, can also provide accurate prediction to experimental values (r² = 0.84). Therefore, the multiple independent sampling method can be more efficient than a single, long simulation method. Since future scaffold expansions may significantly change the benzimidazole's physiochemical properties (charges, etc.) and possibly binding modes, which may affect the sensitivities of various parameters, the relatively insensitive “multiple independent sampling method” may avoid the need of an entirely new validation study. Moreover, due to large fluctuating entropy values, (QM/)MM‐P(G)BSA were limited to inhibitors’ relative affinity prediction, but not the absolute affinity. The developed protocol will support an ongoing benzimidazole lead optimization program. © 2015 Wiley Periodicals, Inc.     

Binding-affinity predictions of HSP90 in the D3R Grand Challenge 2015 with docking,MM/GBSA,QM/MM,and free-energy simulations

We have estimated the binding affinity of three sets of ligands of the heat-shock protein 90 in the D3R grand challenge blind test competition. We have employed four different methods, based on five different crystal structures: first, we docked the ligands to the proteins with induced-fit docking with the Glide software and calculated binding affinities with three energy functions. Second, the docked structures were minimised in a continuum solvent and binding affinities were calculated with the MM/GBSA method (molecular mechanics combined with generalised Born and solvent-accessible surface area solvation). Third, the docked structures were re-optimised by combined quantum mechanics and molecular mechanics (QM/MM) calculations. Then, interaction energies were calculated with quantum mechanical calculations employing 970–1160 atoms in a continuum solvent, combined with energy corrections for dispersion, zero-point energy and entropy, ligand distortion, ligand solvation, and an increase of the basis set to quadruple-zeta quality. Fourth, relative binding affinities were estimated by free-energy simulations, using the multi-state Bennett acceptance-ratio approach. Unfortunately, the results were varying and rather poor, with only one calculation giving a correlation to the experimental affinities larger than 0.7, and with no consistent difference in the quality of the predictions from the various methods. For one set of ligands, the results could be strongly improved (after experimental data were revealed) if it was recognised that one of the ligands displaced one or two water molecules. For the other two sets, the problem is probably that the ligands bind in different modes than in the crystal structures employed or that the conformation of the ligand-binding site or the whole protein changes.     

Fast and accurate predictions of binding free energies using MM‐PBSA and MM‐GBSA

In the drug discovery process, accurate methods of computing the affinity of small molecules with a biological target are strongly needed. This is particularly true for molecular docking and virtual screening methods, which use approximated scoring functions and struggle in estimating binding energies in correlation with experimental values. Among the various methods, MM‐PBSA and MM‐GBSA are emerging as useful and effective approaches. Although these methods are typically applied to large collections of equilibrated structures of protein‐ligand complexes sampled during molecular dynamics in water, the possibility to reliably estimate ligand affinity using a single energy‐minimized structure and implicit solvation models has not been explored in sufficient detail. Herein, we thoroughly investigate this hypothesis by comparing different methods for the generation of protein‐ligand complexes and diverse methods for free energy prediction for their ability to correlate with experimental values. The methods were tested on a series of structurally diverse inhibitors of Plasmodium falciparum DHFR with known binding mode and measured affinities. The results showed that correlations between MM‐PBSA or MM‐GBSA binding free energies with experimental affinities were in most cases excellent. Importantly, we found that correlations obtained with the use of a single protein‐ligand minimized structure and with implicit solvation models were similar to those obtained after averaging over multiple MD snapshots with explicit water molecules, with consequent save of computing time without loss of accuracy. When applied to a virtual screening experiment, such an approach proved to discriminate between true binders and decoy molecules and yielded significantly better enrichment curves. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Computation of relative binding free energy for an inhibitor and its analogs binding with Erk kinase using thermodynamic integration MD simulation

In the present study, we carried out thermodynamic integration molecular dynamics simulation for a pair of analogous inhibitors binding with Erk kinase to investigate how computation performs in reproducing the relative binding free energy. The computation with BCC-AM1 charges for ligands gave ?1.1?kcal/mol, deviated from experimental value of ?2.3?kcal/mol by 1.2?kcal/mol, in good agreement with experimental result. The error of computed value was estimated to be 0.5?kcal/mol. To obtain convergence, switching vdw interaction on and off required approximately 10 times more CPU time than switching charges. Residue-based contributions and hydrogen bonding were analyzed and discussed. Furthermore, subsequent simulation using RESP charge for ligand gave ΔΔG of ?1.6?kcal/mol. The computed results are better than the result of ?5.6?kcal/mol estimated using PBSA method in a previous study. Based on these results, we further carried out computations to predict ΔΔG for five new analogs, focusing on placing polar and nonpolar functional groups at the meta site of benzene ring shown in the Fig.?1, to see if these ligands have better binding affinity than the above ligands. The computations resulted that a ligand with polar –OH group has better binding affinity than the previous examined ligand by ~2.0?kcal/mol and two other ligands have better affinity by ~1.0?kcal/mol. The predicted better inhibitors of this kind should be of interest to experimentalist for future experimental enzyme and/or cell assays.     

FEW: A workflow tool for free energy calculations of ligand binding

In the later stages of drug design projects, accurately predicting relative binding affinities of chemically similar compounds to a biomolecular target is of utmost importance for making decisions based on the ranking of such compounds. So far, the extensive application of binding free energy approaches has been hampered by the complex and time‐consuming setup of such calculations. We introduce the free energy workflow (FEW) tool that facilitates setup and execution of binding free energy calculations with the AMBER suite for multiple ligands. FEW allows performing free energy calculations according to the implicit solvent molecular mechanics (MM‐PB(GB)SA), the linear interaction energy, and the thermodynamic integration approaches. We describe the tool's architecture and functionality and demonstrate in a show case study on Factor Xa inhibitors that the time needed for the preparation and analysis of free energy calculations is considerably reduced with FEW compared to a fully manual procedure. © 2013 Wiley Periodicals, Inc.     

Importance of CH/π hydrogen bonds in recognition of the core motif in proline-recognition domains: an ab initio fragment molecular orbital study

We examined CH/π hydrogen bonds in protein/ligand complexes involving at least one proline residue using the ab initio fragment molecular orbital (FMO) method and the program CHPI. FMO calculations were carried out at the Hartree–Fock (HF)/6‐31G*, HF/6‐31G**, second‐order Møller–Plesset perturbation (MP2)/6‐31G*, and MP2/6‐31G** levels for three Src homology 3 (SH3) domains and five proline‐recognition domains (PRDs) complexed with their corresponding ligand peptides. PRDs use a conserved set of aromatic residues to recognize proline‐rich sequences of specific ligands. Many CH/π hydrogen bonds were identified in these complexes. CH/π hydrogen bonds occurred, in particular, in the central part of the proline‐rich motifs. Our results suggest that CH/π hydrogen bonds are important in the recognition of SH3 and PRDs by their ligand peptides and play a vital role in the signal transduction system. Combined use of the FMO method and CHPI analysis is a valuable tool for the study of protein/protein and protein/ligand interactions and may be useful in rational drug design. © 2011 Wiley Periodicals, Inc. J Comput Chem 2011     

Understanding the basis of I50V‐induced affinity decrease in HIV‐1 protease via molecular dynamics simulations using polarized force field

Human immunodeficiency virus (HIV)‐1 protease is one of the most promising drug target commonly utilized to combat Acquired Immune Deficiency Syndrome (AIDS). However, with the emergence of drug resistance arising from mutations, the efficiency of protease inhibitors (PIs) as a viable treatment for AIDS has been greatly reduced. I50V mutation as one of the most significant mutations occurring in HIV‐1 protease will be investigated in this study. Molecular dynamics (MD) simulation was utilized to examine the effect of I50V mutation on the binding of two PIs namely indinavir and amprenavir to HIV‐1 protease. Prior to the simulations conducted, the electron density distributions of the PI and each residue in HIV‐1 protease are derived by combining quantum fragmentation approach molecular fractionation with conjugate caps and Poisson–Boltzmann solvation model based on polarized protein‐specific charge scheme. The atomic charges of the binding complex are subsequently fitted using delta restrained electrostatic potential (delta‐RESP) method to overcome the poor charge determination of buried atom. This way, both intraprotease polarization and the polarization between protease and the PI are incorporated into partial atomic charges. Through this study, the mutation‐induced affinity variations were calculated and significant agreement between experiments and MD simulations conducted was observed for both HIV‐1 protease‐drug complexes. In addition, the mechanism governing the decrease in the binding affinity of PI in the presence of I50V mutation was also explored to provide insights pertaining to the design of the next generation of anti‐HIV drugs. © 2015 Wiley Periodicals, Inc.     

The F130L mutation in streptavidin reduces its binding affinity to biotin through electronic polarization effect

Recently, Baugh et al. discovered that a distal point mutation (F130L) in streptavidin causes no distinct variation to the structure of the binding pocket but a 1000‐fold reduction in biotin binding affinity. In this work, we carry out molecular dynamics simulations and apply an end‐state free energy method to calculate the binding free energies of biotin to wild type streptavidin and its F130L mutant. The absolute binding affinities based on AMBER charge are repulsive, and the mutation induced binding loss is underestimated. When using the polarized protein‐specific charge, the absolute binding affinities are significantly enhanced. In particular, both the absolute and relative binding affinities are in line with the experimental measurements. Further investigation indicates that polarization effect is indispensable in both the generation of structural ensembles and the calculation of interaction energies. This work verifies Baugh's conjecture that electrostatic polarization effect plays an essential role in modulating the binding affinity of biotin to the streptavidin through F130L mutation. © 2013 Wiley Periodicals, Inc.