Experimental Binding Free Energies Research Articles

Accurate protein-ligand binding free energy estimation using QM/MM on multi-conformers predicted from classical mining minima

Accurate prediction of binding free energy is crucial for the rational design of drug candidates and understanding protein-ligand interactions. To address this, we have developed four protocols that combine QM/MM calculations and the mining minima (M2) method, tested on 9 targets and 203 ligands. Our protocols carry out free energy processing with or without conformational search on the selected conformers obtained from M2 calculations, where their force field atomic charge parameters are substituted with those obtained from a QM/MM calculation. The method achieved a high Pearson’s correlation coefficient (0.81) with experimental binding free energies across diverse targets, demonstrating its generality. Using a differential evolution algorithm with a universal scaling factor of 0.2, we achieved a low mean absolute error of 0.60 kcal mol-1. This performance surpasses many existing methods and is comparable to popular relative binding free energy techniques but at significantly lower computational cost.

Simulating Metal-Imidazole Complexes.

One commonly observed binding motif in metalloproteins involves the interaction between a metal ion and histidine's imidazole side chains. Although previous imidazole-M(II) parameters established the flexibility and reliability of the 12-6-4 Lennard-Jones (LJ)-type nonbonded model by simply tuning the ligating atom's polarizability, they have not been applied to multiple-imidazole complexes. To fill this gap, we systematically simulate multiple-imidazole complexes (ranging from one to six) for five metal ions (Co(II), Cu(II), Mn(II), Ni(II), and Zn(II)) which commonly appear in metalloproteins. Using extensive (40 ns per PMF window) sampling to assemble free energy association profiles (using OPC water and standard HID imidazole charge models from AMBER) and comparing the equilibrium distances to DFT calculations, a new set of parameters was developed to focus on energetic and geometric features of multiple-imidazole complexes. The obtained free energy profiles agree with the experimental binding free energy and DFT calculated distances. To validate our model, we show that we can close the thermodynamic cycle for metal-imidazole complexes with up to six imidazole molecules in the first solvation shell. Given the success in closing the thermodynamic cycles, we then used the same extended sampling method for six other metal ions (Ag(I), Ca(II), Cd(II), Cu(I), Fe(II), and Mg(II)) to obtain new parameters. Since these new parameters can reproduce the one-imidazole geometry and energy accurately, we hypothesize that they will reasonably predict the binding free energy of higher-level coordination numbers. Hence, we did not extend the analysis of these ions up to six imidazole complexes. Overall, the results shed light on metal-protein interactions by emphasizing the importance of ligand-ligand interaction and metal-π-stacking within metalloproteins.

Combining the Fragment Molecular Orbital and GRID Approaches for the Prediction of Ligand-Metalloenzyme Binding Affinity: The Case Study of hCA II Inhibitors.

Polarization and charge-transfer interactions play an important role in ligand-receptor complexes containing metals, and only quantum mechanics methods can adequately describe their contribution to the binding energy. In this work, we selected a set of benzenesulfonamide ligands of human Carbonic Anhydrase II (hCA II)-an important druggable target containing a Zn2+ ion in the active site-as a case study to predict the binding free energy in metalloprotein-ligand complexes and designed specialized computational methods that combine the ab initio fragment molecular orbital (FMO) method and GRID approach. To reproduce the experimental binding free energy in these systems, we adopted a machine-learning approach, here named formula generator (FG), considering different FMO energy terms, the hydrophobic interaction energy (computed by GRID) and logP. The main advantage of the FG approach is that it can find nonlinear relations between the energy terms used to predict the binding free energy, explicitly showing their mathematical relation. This work showed the effectiveness of the FG approach, and therefore, it might represent an important tool for the development of new scoring functions. Indeed, our scoring function showed a high correlation with the experimental binding free energy (R2 = 0.76-0.95, RMSE = 0.34-0.18), revealing a nonlinear relation between energy terms and highlighting the relevant role played by hydrophobic contacts. These results, along with the FMO characterization of ligand-receptor interactions, represent important information to support the design of new and potent hCA II inhibitors.

Robust Prediction of Relative Binding Energies for Protein–Protein Complex Mutations Using Free Energy Perturbation Calculations

Computational free energy-based methods have the potential to significantly improve throughput and decrease costs of protein design efforts. Such methods must reach a high level of reliability, accuracy, and automation to be effectively deployed in practical industrial settings in a way that impacts protein design projects. Here, we present a benchmark study for the calculation of relative changes in protein–protein binding affinity for single point mutations across a variety of systems from the literature, using free energy perturbation (FEP+) calculations. We describe a method for robust treatment of alternate protonation states for titratable amino acids, which yields improved correlation with and reduced error compared to experimental binding free energies. Following careful analysis of the largest outlier cases in our dataset, we assess limitations of the default FEP+ protocols and introduce an automated script which identifies probable outlier cases that may require additional scrutiny and calculates an empirical correction for a subset of charge-related outliers. Through a series of three additional case study systems, we discuss how Protein FEP+ can be applied to real-world protein design projects, and suggest areas of further study.

Efficient and accurate binding free energy calculation of Aβ9-40 protofilament propagation.

Self-assembled aggregation of peptides and proteins into regular amyloid fibrils is associated with several neurodegenerative diseases. In case of Alzheimer's disease proteolytic cleavage products of the amyloid precursor protein form pathological amyloid-beta fibrils in a nucleation and propagation phase. The molecular details and thermodynamic driving forces of amyloid formation are not well understood, but are of high relevance for potential pharmacological interference. We used atomistic binding free energy simulations to calculate the free energy of protofilament propagation by an additional Aβ9-40 peptide binding to the protofilament tip. It requires sampling of relevant conformational transitions which is challenging since the monomeric Aβ9-40 peptide is intrinsically disordered. However, the convergence of umbrella simulations can be enhanced by applying additional restraining potentials on the axial, orientational and conformational degrees of freedom. The improved convergence leads to a much closer agreement with experimental binding free energy data compared to unrestrained umbrella sampling. Moreover, the restraining approach results in a separation of contributions to the total binding free energy. The calculated contributions indicate that the free energy change associated with the restriction of conformational freedom upon propagation makes a large opposing contribution of higher magnitude than the total binding free energy. Finally, optimization of the approach leads to further significant reduction of the computational demand which is crucial for systematic studies on mutations, denaturants and inhibitors in the fibril propagation step.

Computational Workflow for Refining AlphaFold Models in Drug Design Using Kinetic and Thermodynamic Binding Calculations: A Case Study for the Unresolved Inactive Human Adenosine A3 Receptor.

A structure-based drug design pipeline that considers both thermodynamic and kinetic binding data of ligands against a receptor will enable the computational design of improved drug molecules. For unresolved GPCR-ligand complexes, a workflow that can apply both thermodynamic and kinetic binding data in combination with alpha-fold (AF)-derived or other homology models and experimentally resolved binding modes of relevant ligands in GPCR-homologs needs to be tested. Here, as test case, we studied a congeneric set of ligands that bind to a structurally unresolved G protein-coupled receptor (GPCR), the inactive human adenosine A3 receptor (hA3R). We tested three available homology models from which two have been generated from experimental structures of hA1R or hA2AR and one model was a multistate alphafold 2 (AF2)-derived model. We applied alchemical calculations with thermodynamic integration coupled with molecular dynamics (TI/MD) simulations to calculate the experimental relative binding free energies and residence time (τ)-random accelerated MD (τ-RAMD) simulations to calculate the relative residence times (RTs) for antagonists. While the TI/MD calculations produced, for the three homology models, good Pearson correlation coefficients, correspondingly, r = 0.74, 0.62, and 0.67 and mean unsigned error (mue) values of 0.94, 1.31, and 0.81 kcal mol-1, the τ-RAMD method showed r = 0.92 and 0.52 for the first two models but failed to produce accurate results for the multistate AF2-derived model. With subsequent optimization of the AF2-derived model by reorientation of the side chain of R1735.34 located in the extracellular loop 2 (EL2) that blocked ligand's unbinding, the computational model showed r = 0.84 for kinetic data and improved performance for thermodynamic data (r = 0.81, mue = 0.56 kcal mol-1). Overall, after refining the multistate AF2 model with physics-based tools, we were able to show a strong correlation between predicted and experimental ligand relative residence times and affinities, achieving a level of accuracy comparable to an experimental structure. The computational workflow used can be applied to other receptors, helping to rank candidate drugs in a congeneric series and enabling the prioritization of leads with stronger binding affinities and longer residence times.

Thermodynamics of Metal-Acetate Interactions.

Metal ions play crucial roles in protein- and ligand-mediated interactions. They not only act as catalysts to facilitate biological processes but are also important as protein structural elements. Accurately predicting metal ion interactions in computational studies has always been a challenge, and various methods have been suggested to improve these interactions. One such method is the 12-6-4 Lennard-Jones (LJ)-type nonbonded model. Using this model, it has been possible to successfully reproduce the experimental properties of metal ions in aqueous solution. The model includes induced dipole interactions typically ignored in the standard 12-6 LJ nonbonded model. In this we expand the applicability of this model to metal ion-carboxylate interactions. Using 12-6-4 parameters that reproduce the solvation free energies of the metal ions leads to an overestimation of metal ion-acetate interactions, thus, prompting us to fine-tune the model to specifically handle the latter. We also show that the standard 12-6 LJ model significantly falls short in reproducing the experimental binding free energy between acetate and 11 metal ions (Ni(II), Mg(II), Cu(II), Zn(II), Co(II), Cu(I), Fe(II), Mn(II), Cd(II), Ca(II), and Ag(I)). In this study, we describe optimized C4 parameters for the 12-6-4 LJ nonbonded model to be used with three widely employed water models (Transferable Intermolecular Potential with 3 Points (TIP3P), Simple Point Charge Extended (SPC/E), and Optimal Point Charge (OPC) water models). These parameters can accurately match the experimental binding free energy between 11 metal ions and acetate. These parameters can be applied to the study of metalloproteins and transition metal ion channels and transporters, as acetate serves as a representative of the negatively charged amino acid side chains from aspartate and glutamate.

Tuning Potential Functions to Host-Guest Binding Data.

Software to more rapidly and accurately predict protein-ligand binding affinities is of high interest for early-stage drug discovery, and physics-based methods are among the most widely used technologies for this purpose. The accuracy of these methods depends critically on the accuracy of the potential functions that they use. Potential functions are typically trained against a combination of quantum chemical and experimental data. However, although binding affinities are among the most important quantities to predict, experimental binding affinities have not to date been integrated into the experimental data set used to train potential functions. In recent years, the use of host-guest complexes as simple and tractable models of binding thermodynamics has gained popularity due to their small size and simplicity, relative to protein-ligand systems. Host-guest complexes can also avoid ambiguities that arise in protein-ligand systems such as uncertain protonation states. Thus, experimental host-guest binding data are an appealing additional data type to integrate into the experimental data set used to optimize potential functions. Here, we report the extension of the Open Force Field Evaluator framework to enable the systematic calculation of host-guest binding free energies and their gradients with respect to force field parameters, coupled with the curation of 126 host-guest complexes with available experimental binding free energies. As an initial application of this novel infrastructure, we optimized generalized Born (GB) cavity radii for the OBC2 GB implicit solvent model against experimental data for 36 host-guest systems. This refitting led to a dramatic improvement in accuracy for both the training set and a separate test set with 90 additional host-guest systems. The optimized radii also showed encouraging transferability from host-guest systems to 59 protein-ligand systems. However, the new radii are significantly smaller than the baseline radii and lead to excessively favorable hydration free energies (HFEs). Thus, users of the OBC2 GB model currently may choose between GB cavity radii that yield more accurate binding affinities and GB cavity radii that yield more accurate HFEs. We suspect that achieving good accuracy on both will require more far-reaching adjustments to the GB model. We note that binding free-energy calculations using the OBC2 model in OpenMM gain about a 10× speedup relative to corresponding explicit solvent calculations, suggesting a future role for implicit solvent absolute binding free-energy (ABFE) calculations in virtual compound screening. This study proves the principle of using host-guest systems to train potential functions that are transferrable to protein-ligand systems and provides an infrastructure that enables a range of applications.

SARS-CoV-2 omicron RBD forms a weaker binding affinity to hACE2 compared to Delta RBD in in-silico studies

The COVID-19 pandemic sparked an unprecedented race in biotechnology in a search for effective therapies and a preventive vaccine. The continued appearance of SARS-CoV-2 variants of concern (VoCs) further swept the world. The entry of SARS-CoV-2 into cells is mediated by binding the receptor-binding domain (RBD) of the S protein to the cell-surface receptor, human angiotensin-converting enzyme 2 (hACE2). In this study, using a coarse-grained force field to parameterize the system, we employed steered-molecular dynamics (SMD) simulations to reveal the binding of SARS-CoV-2 Delta/Omicron RBD to hACE2. Our benchmarked results demonstrate a good correlation between computed rupture force and experimental binding free energy for known protein-protein systems. Moreover, our findings show that the Omicron RBD has a weaker binding affinity to hACE2, consistent with the respective experimental results. This indicates that our method can effectively be applied to other emerging SARS-CoV-2 strains. Communicated by Ramaswamy H. Sarma

Revisiting MMPBSA by Adoption of MC-Based Surface Area/Volume, ANI-ML Potentials, and Two-Valued Interior Dielectric Constant.

Here, we report the accuracy improvements of molecular mechanics Poisson-Boltzmann surface area (MMPBSA) calculations by adoption of ANI-ML potentials in replacement of MM terms, the use of solvent-accessible surface area (SASA) and volume (SAV) values from the Monte Carlo sampling of the probe, and introducing two different interior dielectric constants for electrostatic interactions of protein-ligand (P-L) and polar solvation term in the MMPBSA calculations. Our results show that the Pearson correlation coefficients of MMPBSA-calculated values with respect to experimental binding free energies can be drastically improved from 0.48 to 0.90 by adoption of ANI-ML potentials in replacement of MM energy terms in the equation, referred to as ANI-PBSA. Moreover, we show that the SASA/SAV-combined equation in the scaled particle theory (SPT) can be a better choice to model nonpolar solvation term, reaching nearly the same accuracy by ANI-PBSA calculations. Finally, we introduce two different values of interior dielectric constants, which could be an alternative strategy between the single and variable constant definitions.

Accelerated Enveloping Distribution Sampling to Probe the Presence of Water Molecules.

Determining the presence of water molecules at protein-ligand interfaces is still a challenging task in free-energy calculations. The inappropriate placement of water molecules results in the stabilization of wrong conformational orientations of the ligand. With classical alchemical perturbation methods, such as thermodynamic integration (TI), it is essential to know the amount of water molecules in the active site of the respective ligands. However, the resolution of the crystal structure and the correct assignment of the electron density do not always lead to a clear placement of water molecules. In this work, we apply the one-step perturbation method named accelerated enveloping distribution sampling (AEDS) to determine the presence of water molecules in the active site by probing them in a fast and straightforward way. Based on these results, we combined the AEDS method with standard TI to calculate accurate binding free energies in the presence of buried water molecules. The main idea is to perturb the water molecules with AEDS such that they are allowed to alternate between regular water molecules and non-interacting dummy particles while treating the ligand with TI over an alchemical pathway. We demonstrate the use of AEDS to probe the presence of water molecules for six different test systems. For one of these, previous calculations showed difficulties to reproduce the experimental binding free energies, and here, we use the combined TI-AEDS approach to tackle these issues.

PDA-Pred: Predicting the binding affinity of protein-DNA complexes using machine learning techniques and structural features.

Protein-DNA interactions play an important role in various biological processes such as gene expression, replication, and transcription. Understanding the important features that dictate the binding affinity of protein-DNA complexes and predicting their affinities is important for elucidating their recognition mechanisms. In this work, we have collected the experimental binding free energy (ΔG) for a set of 391 Protein-DNA complexes and derived several structure-based features such as interaction energy, contact potentials, volume and surface area of binding site residues, base step parameters of the DNA and contacts between different types of atoms. Our analysis on relationship between binding affinity and structural features revealed that the important factors mainly depend on the number of DNA strands as well as functional and structural classes of proteins. Specifically, binding site properties such as number of atom contacts between the DNA and protein, volume of protein binding sites and interaction-based features such as interaction energies and contact potentials are important to understand the binding affinity. Further, we developed multiple regression equations for predicting the binding affinity of protein-DNA complexes belonging to different structural and functional classes. Our method showed an average correlation and mean absolute error of 0.78 and 0.98kcal/mol, respectively, between the experimental and predicted binding affinities on a jack-knife test. We have developed a webserver, PDA-PreD (Protein-DNA Binding affinity predictor), for predicting the affinity of protein-DNA complexes and it is freely available at https://web.iitm.ac.in/bioinfo2/pdapred/.

Fully automated screening of compound libraries in drug discovery using docking and molecular dynamics.

In the initial stage of drug discovery, hundreds of thousands of compounds are screened in high-throughput assays to find a small number of hits. To save time and resources, we develop a computational workflow that ranks drug candidates by their binding affinity to the target protein. First, all compounds in a library are docked to the protein. After this filtering step, the drug candidates with the poses scoring best are grouped into pairs, and their relative free energy of binding to the protein is calculated using the accelerated weight histogram method (AWH) for alchemical transformations. Thanks to random moves between λ-windows, a single AWH simulation samples the complete alchemical path connecting two ligands. The probability to visit a λ-window depends on its estimated free energy such that the sampling of the alchemical transformation and the free energy estimate improve each other in a positive feedback loop. To fully integrate docking and free energy calculations, we implemented routines that generate hybrid topologies for drug candidates without human intervention, i.e., they assign force field parameters (AMBER and CHARMM) to the compound and determine the smallest set of atoms to be alchemically transformed for a given compound pair. Being able to convert docking poses into input files for the molecular dynamics simulations in an automated fashion, we extensively test AWH's robustness to (sub-optimal) initial configurations on a benchmark with 13 proteins and approximately 500 ligand pairs for which experimental relative binding free energies are available. In return, free energies obtained from MD for protein-ligand complexes in the PDBbind data set are used as training data for a machine-learning model improving the scoring of docking poses such that both docking and free energy calculations benefit from integration in the workflow.

Evaluating the ability of end-point methods to predict the binding affinity tendency of protein kinase inhibitors.

Because of the high economic cost of exploring the experimental impact of mutations occurring in kinase proteins, computational approaches have been employed as alternative methods for evaluating the structural and energetic aspects of kinase mutations. Among the main computational methods used to explore the affinity linked to kinase mutations are docking procedures and molecular dynamics (MD) simulations combined with end-point methods or alchemical methods. Although it is known that end-point methods are not able to reproduce experimental binding free energy (ΔG) values, it is also true that they are able to discriminate between a better or a worse ligand through the estimation of ΔG. In this contribution, we selected ten wild-type and mutant cocrystallized EGFR-inhibitor complexes containing experimental binding affinities to evaluate whether MMGBSA or MMPBSA approaches can predict the differences in affinity between the wild type and mutants forming a complex with a similar inhibitor. Our results show that a long MD simulation (the last 50 ns of a 100 ns-long MD simulation) using the MMGBSA method without considering the entropic components reproduced the experimental affinity tendency with a Pearson correlation coefficient of 0.779 and an R2 value of 0.606. On the other hand, the correlation between theoretical and experimental ΔΔG values indicates that the MMGBSA and MMPBSA methods are helpful for obtaining a good correlation using a short rather than a long simulation period.

Computational Analysis of Triazole-Based Kojic Acid Analogs as Tyrosinase Inhibitors by Molecular Dynamics and Free Energy Calculations

Molecular docking, molecular dynamics (MD) simulations and the linear interaction energy (LIE) method were used here to predict binding modes and free energy for a set of 1,2,3-triazole-based KA analogs as potent inhibitors of Tyrosinase (TYR), a key metalloenzyme of the melanogenesis process. Initially, molecular docking calculations satisfactorily predicted the binding mode of evaluated KA analogs, where the KA part overlays the crystal conformation of the KA inhibitor into the catalytic site of TYR. The MD simulations were followed by the LIE method, which reproduced the experimental binding free energies for KA analogs with an r2 equal to 0.97, suggesting the robustness of our theoretical model. Moreover, the van der Waals contributions performed by some residues such as Phe197, Pro201, Arg209, Met215 and Val218 are responsible for the binding recognition of 1,2,3-triazole-based KA analogs in TYR catalytic site. Finally, our calculations provide suitable validation of the combination of molecular docking, MD, and LIE approaches as a powerful tool in the structure-based drug design of new and potent TYR inhibitors.

Dual A1/A3 Adenosine Receptor Antagonists: Binding Kinetics and Structure-Activity Relationship Studies Using Mutagenesis and Alchemical Binding Free Energy Calculations.

Drugs targeting adenosine receptors (AR) can provide treatment for diseases. We report the identification of 7-(phenylamino)-pyrazolo[3,4-c]pyridines L2-L10, A15, and A17 as low-micromolar to low-nanomolar A1R/A3R dual antagonists, with 3-phenyl-5-cyano-7-(trimethoxyphenylamino)-pyrazolo[3,4-c]pyridine (A17) displaying the highest affinity at both receptors with a long residence time of binding, as determined using a NanoBRET-based assay. Two binding orientations of A17 produce stable complexes inside the orthosteric binding area of A1R in molecular dynamics (MD) simulations, and we selected the most plausible orientation based on the agreement with alanine mutagenesis supported by affinity experiments. Interestingly, for drug design purposes, the mutation of L2506.51 to alanine increased the binding affinity of A17 at A1R. We explored the structure-activity relationships against A1R using alchemical binding free energy calculations with the thermodynamic integration coupled with the MD simulation (TI/MD) method, applied on the whole G-protein-coupled receptor-membrane system, which showed a good agreement (r = 0.73) between calculated and experimental relative binding free energies.

Relative binding free energy calculations with transformato: A molecular dynamics engine-independent tool.

We present the software package transformato for the setup of large-scale relative binding free energy calculations. Transformato is written in Python as an open source project (https://github.com/wiederm/transformato); in contrast to comparable tools, it is not closely tied to a particular molecular dynamics engine to carry out the underlying simulations. Instead of alchemically transforming a ligand L 1 directly into another L 2, the two ligands are mutated to a common core. Thus, while dummy atoms are required at intermediate states, in particular at the common core state, none are present at the physical endstates. To validate the method, we calculated 76 relative binding free energy differences for five protein–ligand systems. The overall root mean squared error to experimental binding free energies is 1.17 kcal/mol with a Pearson correlation coefficient of 0.73. For selected cases, we checked that the relative binding free energy differences between pairs of ligands do not depend on the choice of the intermediate common core structure. Additionally, we report results with and without hydrogen mass reweighting. The code currently supports OpenMM, CHARMM, and CHARMM/OpenMM directly. Since the program logic to choose and construct alchemical transformation paths is separated from the generation of input and topology/parameter files, extending transformato to support additional molecular dynamics engines is straightforward.

Accurate Binding Free Energy Method from End-State MD Simulations.

Herein, we introduce a new strategy to estimate binding free energies using end-state molecular dynamics simulation trajectories. The method is adopted from linear interaction energy (LIE) and ANI-2x neural network potentials (machine learning) for the atomic simulation environment (ASE). It predicts the single-point interaction energies between ligand–protein and ligand–solvent pairs at the accuracy of the wb97x/6-31G* level for the conformational space that is sampled by molecular dynamics (MD) simulations. Our results on 54 protein–ligand complexes show that the method can be accurate and have a correlation of R = 0.87–0.88 to the experimental binding free energies, outperforming current end-state methods with reduced computational cost. The method also allows us to compare BFEs of ligands with different scaffolds. The code is available free of charge (documentation and test files) at https://github.com/otayfuroglu/deepQM.

Hazardous Shortcuts in Standard Binding Free Energy Calculations.

Calculating the standard binding free energies of protein-protein and protein-ligand complexes from atomistic molecular dynamics simulations in explicit solvent is a problem of central importance in computational biophysics. A rigorous strategy for carrying out such calculations is the so-called "geometrical route". In this method, two molecular objects are progressively separated from one another in the presence of orientational and conformational restraints serving to control the change in configurational entropy that accompanies the dissociation process, thereby allowing the computations to converge within simulations of affordable length. Although the geometrical route provides a rigorous theoretical framework, a tantalizing computational shortcut consists of simply leaving out such orientational and conformational degrees of freedom during the separation process. Here the accuracy and convergence of the two approaches are critically compared in the case of two protein-ligand complexes (Abl kinase-SH3:p41 and MDM2-p53:NVP-CGM097) and three protein-protein complexes (pig insulin dimer, SARS-CoV-2 spike RBD:ACE2, and CheA kinase-P2:CheY). The results of the simulations that strictly follow the geometrical route match the experimental standard binding free energies within chemical accuracy. In contrast, simulations bereft of geometrical restraints converge more poorly, yielding inconsistent results that are at variance with the experimental measurements. Furthermore, the orientational and positional time correlation functions of the protein in the unrestrained simulations decay over several microseconds, a time scale that is far longer than the typical simulation times of the geometrical route, which explains why those simulations fail to sample the relevant degrees of freedom during the separation process of the complexes.

Protein-ligand binding affinity prediction of cyclin-dependent kinase-2 inhibitors by dynamically averaged fragment molecular orbital-based interaction energy.

Fragment molecular orbital (FMO) method is a powerful computational tool for structure-based drug design, in which protein-ligand interactions can be described by the inter-fragment interaction energy (IFIE) and its pair interaction energy decomposition analysis (PIEDA). Here, we introduced a dynamically averaged (DA) FMO-based approach in which molecular dynamics simulations were used to generate multiple protein-ligand complex structures for FMO calculations. To assess this approach, we examined the correlation between the experimental binding free energies and DA-IFIEs of six CDK2 inhibitors whose net charges are zero. The correlation between the experimental binding free energies and snapshot IFIEs for X-ray crystal structures was R2 =0.75. Using the DA-IFIEs, the correlation significantly improved to 0.99. When an additional CDK2 inhibitor with net charge of -1 was added, the DA FMO-based scheme with the dispersion energies still achieved R2 =0.99, whereas R2 decreased to 0.32 employing all the energy terms of PIEDA.