Free Energy Methods Research Articles

Combining QM/MM Calculations with Classical Mining Minima to Predict Protein-Ligand Binding Free Energy.

We developed an effective binding free energy prediction protocol which incorporates quantum mechanical/molecular mechanical (QM/MM) calculations to substitute the specified atomic charges of force fields with quantum-mechanically recalculated ones at a proposed pose using a mining minima approach with the VeraChem mining minima engine. We tested this protocol using seven well-known targets with 147 different ligands and compared it with classical mining minima and the most popular binding free energy (BFE) methods using different metrics. Our new protocol, dubbed Qcharge-VM2, yielded an overall Pearson correlation of 0.86, which was better than all the methods examined. Qcharge-VM2 performed significantly better than implicit solvent-based methods, such as MM-GBSA and MM-PBSA, but not as good as explicit water-based free energy perturbation methods, such as FEP+, in terms of root-mean-square error, RMSE (1.75 kcal/mol) and mean unsigned error, MUE (1.39 kcal/mol) on a limited set of targets. However, our protocol is substantially less computationally demanding compared with FEP+. The combined accuracy and efficiency of our method can be valuable in drug discovery campaigns.

Unveiling MurE ligase potential inhibitors for treating multi-drug resistant Acinetobacter baumannii

Acinetobacter baumannii is an opportunistic pathogen with ability to cause serious infection such as bacteremia, ventilator associated pneumonia, and wound infections. As strains of A. baumannii are resistant to almost all clinically used antibiotics and with the emergence of carbapenems resistant phenotypes warrants the search for novel antibiotics. Considering this, herein, a series of computer aided drug designing approach was utilized to search novel chemical scaffolds that bind stronger to MurE ligase enzyme of A. baumannii, which is involved peptidoglycan synthesis. The work identified LAS_22461675, LAS_34000090 and LAS_51177972 compounds as promising binding molecules with MurE enzyme having binding energy score of −10.5 kcal/mol, −9.3 kcal/mol and −8.6 kcal/mol, respectively. The compounds were found to achieve docked inside the MurE substrate binding pocket and established close distance chemical interactions. The interaction energies were dominated by van der Waals and less contributions were seen from hydrogen bonding energy. The dynamic simulation assay predicted the complexes stable with no major global and local changes noticed. The docked stability was also validated by MM/PBSA and MM/GBSA binding free energy methods. The net MM/GBSA binding free energy of LAS_22461675 complex, LAS_34000090 complex and LAS_51177972 complex is −26.25 kcal/mol, −27.23 kcal/mol and −29.64 kcal/mol, respectively. Similarly in case of MM-PBSA, the net energy value was in following order; LAS_22461675 complex (−27.67 kcal/mol), LAS_34000090 complex (−29.94 kcal/mol) and LAS_51177972 complex (−27.32 kcal/mol). The AMBER entropy and WaterSwap methods also confirmed stable complexes formation. Further, molecular features of the compounds were determined that predicted compounds to have good druglike properties and pharmacokinetic favorable. The study concluded the compounds to good candidates to be tested by in vivo and in vitro experimental assays. Communicated by Ramaswamy H. Sarma

QSAR via Multisite λ-Dynamics in the Orphaned TSSK1B Kinase.

Multisite λ-dynamics (MSλD) is a novel method for the calculation of relative free energies of binding for ligands to their targeted receptors. It can be readily used to examine a large number of molecules with multiple functional groups at multiple sites around a common core. This makes MSλD a powerful tool in structure-based drug design. In the present study, MSλD is applied to calculate the relative binding free energies of 1296 inhibitors to the testis specific serine kinase 1B (TSSK1B), a validated target for male contraception. For this system, MSλD requires significantly fewer computational resources compared to traditional free energy methods like free energy perturbation or thermodynamic integration. From MSλD simulations, we examined whether modifications of a ligand at two different sites are coupled or not. Based on our calculations, we established a quantitative structure-activity relationship (QSAR) for this set of molecules and identified a site in the ligand where further modification, such as adding more polar groups, may lead to increased binding affinity. This article is protected by copyright. All rights reserved.

Accounting for Solvation Correlation Effects on the Thermodynamics of Water Networks in Protein Cavities

Macromolecular recognition and ligand binding are at the core of biological function and drug discovery efforts. Water molecules play a significant role in mediating the protein-ligand interaction, acting as more than just the surrounding medium by affecting the thermodynamics and thus the outcome of the binding process. As individual water contributions are impossible to measure experimentally, a range of computational methods have emerged to identify hydration sites in protein pockets and characterize their energetic contributions for drug discovery applications. Even though several methods model solvation effects explicitly, they focus on determining the stability of specific water sites independently and neglect solvation correlation effects upon replacement of clusters of water molecules, which typically happens in hit-to-lead optimization. In this work, we rigorously determine the conjoint effects of replacing all combinations of water molecules in protein binding pockets through the use of the RE-EDS multistate free-energy method, which combines Hamiltonian replica exchange (RE) and enveloping distribution sampling (EDS). Applications on the small bovine pancreatic trypsin inhibitor and four proteins of the bromodomain family illustrate the extent of solvation correlation effects on water thermodynamics, with the favorability of replacement of the water sites by pharmacophore probes highly dependent on the composition of the water network and the pocket environment. Given the ubiquity of water networks in biologically relevant protein targets, we believe our approach can be helpful for computer-aided drug discovery by providing a pocket-specific and a priori systematic consideration of solvation effects on ligand binding and selectivity.

Elastic properties of SiO2 nanostructure in high-pressure conditions

In this study, the elastic properties of two high-pressure polymorphs SiO2 nanostructure, stishovite and CaCl2-type, are obtained using Density Functional Theory in 0-80 GPa high pressure domain at zero temperature, based on reducing an interacting many-electron problem to a single-electron problem. It is shown that below 40 GPa, the stishovite phase is more stable; superior to this limit, the CaCl2-type phase becomes more stable, using Gibbs free energy method. Furthermore, the pressure dependence of the density, volume, bulk, and shear moduli were defined in the selected pressure domain.

Development and test of highly accurate endpoint free energy methods. 2: Prediction of logarithm of n-octanol-water partition coefficient (logP) for druglike molecules using MM-PBSA method.

The logarithm of n-octanol-water partition coefficient (logP) is frequently used as an indicator of lipophilicity in drug discovery, which has substantial impacts on the absorption, distribution, metabolism, excretion, and toxicity of a drug candidate. Considering that the experimental measurement of the property is costly and time-consuming, it is of great importance to develop reliable prediction models for logP. In this study, we developed a transfer free energy-based logP prediction model-FElogP. FElogP is based on the simple principle that logP is determined by the free energy change of transferring a molecule from water to n-octanol. The underlying physical method to calculate transfer free energy is the molecular mechanics-Poisson Boltzmann surface area (MM-PBSA), thus this method is named as free energy-based logP (FElogP). The superiority of FElogP model was validated by a large set of 707 structurally diverse molecules in the ZINC database for which the measurement was of high quality. Encouragingly, FElogP outperformed several commonly-used QSPR or machine learning-based logP models, as well as some continuum solvation model-based methods. The root-mean-square error (RMSE) and Pearson correlation coefficient(R) between the predicted and measured values are 0.91 log units and 0.71, respectively, while the runner-up, the logP model implemented in OpenBabel had an RMSE of 1.13 log units and R of 0.67. Given the fact that FElogP was not parameterized against experimental logP directly, its excellent performance is likely to be expanded to arbitrary organic molecules covered by the general AMBER force fields.

Development and test of highly accurate endpoint free energy methods. 1: Evaluation of ABCG2 charge model on solvation free energy prediction and optimization of atom radii suitable for more accurate solvationfreeenergy prediction by the PBSA method.

Accurate estimation of solvation free energy (SFE) lays the foundation for accurate prediction of binding free energy. The Poisson-Boltzmann (PB) or generalized Born (GB) combined with surface area (SA) continuum solvation method (PBSA and GBSA) have been widely used in SFE calculations because they can achieve good balance between accuracy and efficiency. However, the accuracy of these methods can be affected by several factors such as the charge models, polar and nonpolar SFE calculation methods and the atom radii used in the calculation. In this work, the performance of the ABCG2 (AM1-BCC-GAFF2) charge model as well as other two charge models, that is, RESP (Restrained Electrostatic Potential) and AM1-BCC (Austin Model 1-bond charge corrections), on the SFE prediction of 544 small molecules in water by PBSA/GBSA was evaluated. In order to improve the performance of the PBSA prediction based on the ABCG2 charge, we further explored the influence of atom radii on the prediction accuracy and yielded a set of atom radius parameters for more accurate SFE prediction using PBSA based on the ABCG2/GAFF2 by reproducing the thermodynamic integration (TI) calculation results. The PB radius parameters of carbon, oxygen, sulfur, phosphorus, chloride, bromide and iodine, were adjusted. New atom types, on, oi, hn1, hn2, hn3, were introduced to further improve the fitting performance. Then, we tuned the parameters in the nonpolar SFE model using the experimental SFE data and the PB calculation results. By adopting the new radius parameters and new nonpolar SFE model, the root mean square error (RMSE) of the SFE calculation for the 544 molecules decreased from 2.38 to 1.05 kcal/mol. Finally, the new radius parameters were applied in the prediction of protein-ligand binding free energies using the MM-PBSA method. For the eight systems tested, we could observe higher correlation between the experiment data and calculation results and smaller prediction errors for the absolute binding free energies, demonstrating that our new radius parameters can improve the free energy calculation using the MM-PBSA method.

The electronic structure of genome editors from the first principles

Ab-initio molecular dynamics enables following the dynamics of biological systems from the first principles, describing the electronic structure and offering the opportunity to “watch” the evolution of biochemical processes with unique resolution, beyond the capabilities of state-of-the-art experimental techniques. This article reports the role of first-principles (ab-initio) molecular dynamics (MD) in the CRISPR-Cas9 genome editing revolution, achieving a profound understanding of the enzymatic function and offering valuable insights for enzyme engineering. We introduce the methodologies and explain the use of ab-initio MD simulations to establish the two-metal dependent mechanism of DNA cleavage in the RuvC domain of the Cas9 enzyme, and how a second catalytic domain, HNH, cleaves the target DNA with the aid of a single metal ion. A detailed description of how ab-initio MD is combined with free-energy methods—i.e., thermodynamic integration and metadynamics—to break and form chemical bonds is given, explaining the use of these methods to determine the chemical landscape and establish the catalytic mechanism in CRISPR-Cas9. The critical role of classical methods is also discussed, explaining theory and application of constant pH MD simulations, used to accurately predict the catalytic residues’ protonation states. Overall, first-principles methods are shown to unravel the electronic structure and reveal the catalytic mechanism of the Cas9 enzyme, providing valuable insights that can serve for the design of genome editing tools with improved catalytic efficiency or controllable activity.

High-throughput binding free energy simulations: Applications in drug discovery.

Computational drug discovery methods are becoming more and more prevalent in the drug discovery pipeline. Binding free energy (BFE) simulations are a class of molecular dynamics-based free energy methods that utilize a nonphysical “alchemical” thermodynamic pathway to efficiently and rigorously compute protein-ligand complex binding free energies. With the advent of GPU computing, BFE simulations can be performed on a timescale reasonable for small molecule drug discovery applications. Combining these methods with machine learning methods such as active learning, one can rapidly explore large chemical spaces with a fraction of the cost of brute force simulations. Herein, we describe a framework for the high-throughput and highly automated screening of a large number of chemically diverse small molecule ligands for the LRRK2 WDR domain using BFE calculations with thermodynamic integration. LRRK2 is the most commonly mutated gene in familial Parkinson's Disease, and its WDR domain is an understudied drug target with no known molecular inhibitors. This work was done as a submission to the Critical Assessment of Computational Hit-Finding Experiments (CACHE) Challenge. We also discuss various methods to minimize computational cost while maintaining accuracy and make recommendations for future high-throughput screening campaigns.

Investigating mutant hydrophobic patch exposure in troponin C using molecular dynamics-based free energy methods.

Contraction of cardiomyocytes is orchestrated through calcium-binding to the N-terminal regulatory domain of cardiac troponin C (NcTnC). This binding event initiates a series of dynamic and conformational changes that expose a hydrophobic patch to which troponin I (TnI) associates. Cardiomyopathies have been linked to mutations in NcTnC that disrupt these precise interactions and modulate free energy of opening. Mutations that lower the free energy of opening have a calcium-sensitizing effect and those that increase free energy of opening have a calcium-desensitizing effect, like mutations implicated with hypertrophic (HCM) and dilated cardiomyopathy (DCM), respectively. Albeit mutations occurring in the loop region of NcTnC have been most extensively studied, information regarding mutations within the hydrophobic patch is lacking. Simulating the closed-to-open transition in NcTnC is challenging and typically requires microsecond timescales just to observe partial opening of the patch. Here, we use targeted molecular dynamics (TMD) to efficiently sample the landscape between the open and closed conformations. Umbrella sampling (US) was employed to predict the free energy of opening by using frames sourced from TMD, and preliminary data corroborates the expected trends for known HCM and DCM mutations in NcTnC. RosettaDesign was used to prepare and score novel mutation candidates by systematically mutating residues at eight different positions in the patch to almost every other amino acid. Based on the ref15_cart scoring function, up to three mutation candidates from each position were selected to proceed to further evaluation with US. In this work, we identify several mutations that are expected to perturb patch opening and present an innovative approach to quantitatively study the influence hydrophobic patch residues play in the molecular basis of myocyte contraction.

Assessing the effect of short and long-term aging on moisture damage of hot mix asphalt using two different methods

ABSTRACT The effects of short-term aging (STA) and long-term aging (LTA) on moisture damage mechanisms were investigated using mechanical and thermodynamic methods. PG58-22 and PG64-16 asphalt binders modified with Dodecylamine (DA) and Wetfix (WT) additives and limestone and granite aggregates were used. The moisture sensitivity was assessed using the indirect tensile strength (ITS) test via the modified Lottman test and the surface free energy (SFE) method. As asphalt binder becomes harder due to aging, ITS values increase in dry conditions. However, ITS values decrease due to the texture and mineralogical properties of aggregates, the reduction of mastic cohesion and asphalt binder-aggregate adhesion because of water infiltration and weakening of the mixture in wet conditions. The tensile strength ratio (TSR) test revealed that TSR values decrease due to aging. The cohesion free energy of the PG64-16 asphalt binder increased by 8.9% in the STA state and by 29.8% in the LTA state compared to the unaged state. This rise in PG58-22 was about the same as that in PG64-16. Asphalt binder aging in all the samples reduced the adhesion energy. These variations indicated that aging reduces the probability of cohesion failure and increases the probability of adhesion failure.

Free Energy Methods for the Description of Molecular Processes.

Efforts to combine theory and experiment to advance our knowledge of molecular processes relevant to biophysics have been considerably enhanced by the contribution of statistical-mechanics simulations. Key to the understanding of such molecular processes is the underlying free-energy change. Being able to accurately predict this change from first principles represents an appealing prospect. Over the past decades, the synergy between steadily growing computational resources and unrelenting methodological developments has brought free-energy calculations into the arsenal of tools commonly utilized to tackle important questions that experiment alone has left unresolved. The continued emergence of new options to determine free energies has also bred confusion amid the community of users, who may find it difficult to choose the best-suited algorithm to address the problem at hand. In an attempt to clarify the current landscape, this review recounts how the field has been shaped and how the broad gamut of methods available today is rooted in a few foundational principles laid down many years ago.Three examples of molecular processes central to biophysics illustrate where free-energy calculations stand and what are the conceptual and practical obstacles that we must overcome to increase their predictive power.

Mechanistic Insights into the Reactive Uptake of Chlorine Nitrate at the Air-Water Interface.

It is well-known that the aqueous-phase processing of chlorine nitrate (ClONO2) plays a crucial role in ozone depletion. However, many of the physical and chemical properties of ClONO2 at the air-water interface or in bulk water are unknown or not understood on a microscopic scale. Here, the solvation and hydrolysis of ClONO2 at the air-water interface and in bulk water at 300 K were investigated by classical and ab initio molecular dynamics (AIMD) simulations combined with free energy methods. Our results revealed that ClONO2 prefers to accumulate at the air-water interface rather than in the bulk phase. Specifically, halogen bonding interactions (ClONO2)Cl···O(H2O) were found to be the predominant interactions between ClONO2 and H2O. Moreover, metadynamics-biased AIMD simulations revealed that ClONO2 hydrolysis is catalyzed at the air-water interface with an activation barrier of only ∼0.2 kcal/mol; additionally, the difference in free energy between the product and reactant is only ∼0.1 kcal/mol. Surprisingly, the near-barrierless reaction and the comparable free energies of the reactant and product suggested that the ClONO2 hydrolysis at the air-water interface is reversible. When the temperature is lowered from 300 to 200 K, the activation barrier for the ClONO2 hydrolysis at the air-water interface is increased to ∼5.4 kcal/mol. These findings have important implications for the interpretation of experiments.

Free energy calculations and unbiased molecular dynamics targeting the liquid-liquid transition in water no man's land.

The existence of a first-order phase transition between a low-density liquid (LDL) and a high-density liquid (HDL) form of supercooled water has been a central and highly debated issue of physics and chemistry for the last three decades. We present a computational study that allows us to determine the free-energy landscapes of supercooled water over a wide range of pressure and temperature conditions using the TIP4P/2005 force field. Our approach combines topology-based structural transformation coordinates, state-of-the-art free-energy calculation methods, and extensive unbiased molecular dynamics. All our diverse simulations cannot detect any barrier within the investigated timescales and system size, for a discontinuous transition between the LDL and HDL forms throughout the so-called "no man's land," until the onset of the solid, non-diffusive amorphous forms.

Uncovering water effects in protein-ligand recognition: importance in the second hydration shell and binding kinetics.

Developing a ligand with high affinity for a specific protein target is essential for drug design, and water molecules are well known to play a key role in protein-drug recognition. However, predicting the role of particularly ordered water molecules in drug binding remains challenging. Furthermore, hydration free energy contributed from the water network, including the second shell of water molecules, is far from being well studied. In this research we focused on these aspects to accurately and efficiently evaluate water effects in protein-ligand binding affinity. We developed a new strategy using a free-energy calculation method, VM2. We successfully predicted the stable ordered water molecules in a number of protein systems: PDE 10a, HSP90, tryptophan synthase (TRPS), CDK2 and Factor Xa. In some of these, the second shell of water molecules appeared to be critical in protein-ligand binding. We also applied the strategy to largely improve binding free energy calculation using the MM/PBSA method. When applying MM/PBSA alone for two systems, CDK2 and Factor Xa, the computed binding free energy resulted in poor to moderate R2 values with experimental data. However, including water free energy correction greatly improved the free energy calculation. Furthermore, our work helped to explain how xk263 is a 1000 times faster binder to HIVp than ritonavir, a potentially useful tool for investigating binding kinetics. Our studies reveal the importance of fully considering water effects in therapeutic developments in pharmaceutical and biotechnology industries and for fundamental research in protein-ligand recognition.

Solvation free energies of alcohols in water: temperature and pressure dependences.

Solvation free energies μ* of amphiphilic species, methanol and 1,2-hexanediol, are obtained as a function of temperature or pressure based on molecular dynamics simulations combined with efficient free-energy calculation methods. In general, μ* of an amphiphile can be divided into and , the nonpolar and electrostatic contributions, and the former is further divided into and which are the work of cavity formation process and the free energy change due to weak, attractive interactions between the solute molecule and surrounding solvent molecules. We demonstrate that μ* of the two amphiphilic solutes can be obtained accurately using a perturbation combining method, which relies on the exact expressions for and and requires no simulations of intermediate systems between the solute with strong, repulsive interactions and the solute with the van der Waals and electrostatic interactions. The decomposition of μ* gives us several physical insights including that μ* is an increasing function of T due to , that the contributions of hydrophilic groups to the temperature dependence of μ* are additive, and that the contribution of the van der Waals attraction to the solvation volume is greater than that of the electrostatic interactions.

Transport Coefficients of H2–Cu/Al2O3 Plasma Mixtures Considering Solid and Liquid Status

Obtaining the plasma transport coefficients of the mixed system of arc-extinguishing medium and erosion product is the precondition of numerical simulation of erosion–diffusion in an arc chamber. In this article, the minimum Gibbs free energy method was used to calculate the composition of the H2–Cu/Al2O3 system, including solid and liquid copper (aluminum oxide), because copper and aluminum oxide are commonly used as electrode and ceramic materials in hydrogen dc contactors, respectively. The transport coefficients of H2–Cu/Al2O3 system and the phase transition of erosion products provide the basis for the subsequent simulation of erosion–diffusion of electrode and arc chamber in the H2 plasma. Results revealed that: 1) in the low-temperature range, copper (aluminum oxide) mainly exists in the form of a solid or a liquid state, resulting in the transport coefficients of the H2–Cu/Al2O3 system consistent with that of the H2 system; 2) as temperature increases, the copper (aluminum oxide) undergoes gasification (sublimation–decomposition) reaction that exists in a gaseous state and begins to affect the transport properties of the H2–Cu/Al2O3 system; in particular, the proportion of copper (aluminum oxide) has the most significant effect on viscosity; and 3) as the temperature continues to rise, atoms and molecules in the H2–Cu/Al2O3 system are completely ionized and the proportion of copper (aluminum oxide) has little effect on the transport properties (except electrical conductivity) of the H2–Cu/Al2O3 system. Finally, the transport coefficients of H2–Cu/Al2O3 plasma mixtures are reliable enough, in engineering application aspect, to be applied for simulation of erosion–diffusion in hydrogen dc contactors.

Accelerated Enveloping Distribution Sampling (AEDS) Allows for Efficient Sampling of Orthogonal Degrees of Freedom.

One of the most challenging aspects in the molecular simulation of proteins is the study of slowly relaxing processes, such as loop rearrangements or ligands that adopt different conformations in the binding site. State-of-the-art methods used to calculate binding free energies rely on performing several short simulations (lambda steps), in which the ligand is slowly transformed into the endstates of interest. This makes capturing the slowly relaxing processes even more difficult, as they would need to be observed in all of the lambda steps. One attractive alternative is the use of a reference state capable of sampling all of the endstates of interest in a single simulation. However, the energy barriers between the states can be too high to overcome, thus hindering the sampling of all of the relevant conformations. Accelerated enveloping distribution sampling (AEDS) is a recently developed reference state technique that circumvents the high-energy-barrier challenge by adding a boosting potential that flattens the energy landscape without distorting the energy minima. In the present work, we apply AEDS to the well-studied benchmark system T4 lysozyme L99A. The T4 lysozyme L99A mutant contains a hydrophobic pocket in which there is a valine (valine 111), whose conformation influences the binding efficiencies of the different ligands. Incorrectly sampling the dihedral angle can lead to errors in predicted binding free energies of up to 16 kJ mol-1. This protein constitutes an ideal scenario to showcase the advantages and challenges when using AEDS in the presence of slow relaxing processes. We show that AEDS is able to successfully sample the relevant degrees of freedom, providing accurate binding free energies, without the need of previous information of the system in the form of collective variables. Additionally, we showcase the capabilities of AEDS to efficiently screen a set of ligands. These results represent a promising first step toward the development of free-energy methods that can respond to more intricate molecular events.

The Thermal Stability of Asymmetric Separated Configurations inside Alloy Nanoparticles: Atomic-Scale Modeling of Pd-Ir Nanophase Diagrams.

Compared to alloy bulk phase diagrams, the experimental determination of phase diagrams for alloy nanoparticles (NPs), which are useful in various nanotechnological applications, involves significant technical difficulties, making theoretical modeling a feasible alternative. Yet, being quite challenging, modeling of separation nanophase diagrams is scarce in the literature. The task of predicting comprehensive nanophase diagrams for Pd-Ir face-centered cubic-based three cuboctahedra is facilitated in this study by combining the computationally efficient statistical-mechanical Free-energy Concentration Expansion Method, which includes short-range order (SRO) with coordination-dependent bond-energy variations as part of the input and with rotationally symmetric site grouping for extra efficiency. This nanosystem has been chosen mainly because of the very small atomic mismatch that simplifies the modeling, e.g., in the assessment of vibrational entropy contributions based in this work on fitting to the Pd-Ir experimental bulk critical temperature. This entropic effect, together with SRO, leads to significant destabilization of low-T Quasi-Janus (QJ) asymmetric configurations of the NP core, which transform to symmetric partially mixed nanophases. First-order and second-order intracore transitions are predicted for dilute and intermediate-range compositions, respectively. Caloric curves computed for the former case yield the NP-size dependent transition latent heat, and in the latter case critical temperatures exhibit a specific scaling behavior. The computed separation diagrams and intracore solubility diagrams reflect enhanced elemental mixing in smaller QJ nanophases. In addition to these diagrams, the revealed near-surface compositional variations are likely to be pertinent to the utilization of Pd-Ir NPs, e.g., in catalysis.

Prioritizing Small Sets of Molecules for Synthesis through in-silico Tools: A Comparison of Common Ranking Methods.

Prioritizing molecules for synthesis is a key role of computational methods within medicinal chemistry. Multiple tools exist for ranking molecules, from the cheap and popular molecular docking methods to more computationally expensive molecular-dynamics (MD)-based methods. It is often questioned whether the accuracy of the more rigorous methods justifies the higher computational cost and associated calculation time. Here, we compared the performance on ranking the binding of small molecules for seven scoring functions from five docking programs, one end-point method (MM/GBSA), and two MD-based free energy methods (PMX, FEP+). We investigated 16 pharmaceutically relevant targets with a total of 423 known binders. The performance of docking methods for ligand ranking was strongly system dependent. We observed that MD-based methods predominantly outperformed docking algorithms and MM/GBSA calculations. Based on our results, we recommend the application of MD-based free energy methods for prioritization of molecules for synthesis in lead optimization, whenever feasible.