Ligand Binding Free Energy Research Articles

Memprot.GPCR-ModSim: modelling and simulation of membrane proteins in a nutshell.

Memprot.GPCR-ModSim leverages our previous web-based protocol, which was limited to class-A G protein-coupled receptors, to become the first one-stop web server for the modelling and simulation of any membrane protein system. Motivated by the exponential growth of experimental structures and the breakthrough of deep-learning-based structural modelling, the server accepts as input either a membrane-protein sequence, in which case it reports the associated AlphaFold model, or a 3D (experimental, modelled) structure, including quaternary complexes with associated proteins and/or ligands of any kind. In both cases, the molecular dynamics (MD) protocol produces a membrane-embedded, solvated, and equilibrated system, ready to be used as a starting point for further MD simulations, including ligand-binding free energy calculations. Memprot.GPCR-ModSim web server is publicly available at https://memprot.gpcr-modsim.org/. The standalone modules for 3D modelling (PyModSim) or membrane embedding and MD equilibration (PyMemDyn) are available under CC BY-NC 4.0 license terms at the GitHub repository https://github.com/GPCR-ModSim/.

Exploring Binding Pockets in the Conformational States of the SARS-CoV-2 Spike Trimers for the Screening of Allosteric Inhibitors Using Molecular Simulations and Ensemble-Based Ligand Docking.

Understanding mechanisms of allosteric regulation remains elusive for the SARS-CoV-2 spike protein, despite the increasing interest and effort in discovering allosteric inhibitors of the viral activity and interactions with the host receptor ACE2. The challenges of discovering allosteric modulators of the SARS-CoV-2 spike proteins are associated with the diversity of cryptic allosteric sites and complex molecular mechanisms that can be employed by allosteric ligands, including the alteration of the conformational equilibrium of spike protein and preferential stabilization of specific functional states. In the current study, we combine conformational dynamics analysis of distinct forms of the full-length spike protein trimers and machine-learning-based binding pocket detection with the ensemble-based ligand docking and binding free energy analysis to characterize the potential allosteric binding sites and determine structural and energetic determinants of allosteric inhibition for a series of experimentally validated allosteric molecules. The results demonstrate a good agreement between computational and experimental binding affinities, providing support to the predicted binding modes and suggesting key interactions formed by the allosteric ligands to elicit the experimentally observed inhibition. We establish structural and energetic determinants of allosteric binding for the experimentally known allosteric molecules, indicating a potential mechanism of allosteric modulation by targeting the hinges of the inter-protomer movements and blocking conformational changes between the closed and open spike trimer forms. The results of this study demonstrate that combining ensemble-based ligand docking with conformational states of spike protein and rigorous binding energy analysis enables robust characterization of the ligand binding modes, the identification of allosteric binding hotspots, and the prediction of binding affinities for validated allosteric modulators, which is consistent with the experimental data. This study suggested that the conformational adaptability of the protein allosteric sites and the diversity of ligand bound conformations are both in play to enable efficient targeting of allosteric binding sites and interfere with the conformational changes.

Studying the Collective Functional Response of a Receptor in Alchemical Ligand Binding Free Energy Simulations with Accelerated Solvation Layer Dynamics.

Ligand binding free energy simulations (LB-FES) that involve sampling of protein functional conformations have been longstanding challenges in research on molecular recognition. Particularly, modeling of the conformational transition pathway and design of the heuristic biasing mechanism are severe bottlenecks for the existing enhanced configurational sampling (ECS) methods. Inspired by the key role of hydration in regulating conformational dynamics of macromolecules, this report proposes a novel ECS approach that facilitates binding-associated structural dynamics by accelerated hydration transitions in combination with the λ-exchange of free energy perturbation (FEP). Two challenging protein-ligand binding processes involving large configurational transitions of the receptor are studied, with hydration transitions at binding sites accelerated by Hamiltonian-simulated annealing of the hydration layer. Without the need for pathway analysis or ad hoc barrier flattening potential, LB-FES were performed with FEP/λ-exchange molecular dynamics simulation at a minor overhead for annealing of the hydration layer. The LB-FES studies showed that the accelerated rehydration significantly enhances the collective conformational transitions of the receptor, and convergence of binding affinity calculations is obtained at a sweet-spot simulation time scale. Alchemical LB-FES with the proposed ECS strategy is free from the effort of trial and error for the setup and realizes efficient on-the-fly sampling for the collective functional response of the receptor and bound water and therefore presents a practical approach to high-throughput screening in drug discovery.

Precise cuts for tailoring chromene-phenyl COX inhibitors with Ligand Designer

Cyclooxygenases 1 and 2 (COX-1/2) are enzymes renowned for inducing inflammatory responses through the production of prostaglandins. Thus, the development of COX inhibitors has been a promising approach for identifying compounds with anti-inflammatory potential. In this study, we designed 27 new compounds (1–27) based on the structure of a previously known COX inhibitor, using the Ligand Designer tool. Our aim was to improve the affinity of the compounds with COX enzymes by inducing interactions with residue Arg120 while retaining the good π-π stacking interactions of the chromene-phenyl scaffold. Through screening based on ligand-binding free energy defined by molecular docking simulations and MM/GBSA technique, compounds 9 and 10 were identified as having the highest ability to inhibit COX proteins. The binding affinities of the two compounds with COX-1/2 were superior to those of the original NAI10 compound and the reference drug indomethacin. Our virtual screening suggests that compounds 9 and 10 have a strong ability to inhibit COX-1/2 and thus could be promising candidates for further anti-inflammatory drug studies. In essence, our study underscores the pivotal role of the N-aryl iminocoumarin scaffold in shaping the future landscape of novel anti-inflammatory drug development.

Natural compounds inhibit Monkeypox virus methyltransferase VP39 in silico studies

VP39, an essential 2′-O-RNA methyltransferase enzyme discovered in Monkeypox virus (MPXV), plays a vital role in viral RNA replication and transcription. Inhibition of the enzyme may prevent viral replication. In this context, using a combination of molecular docking and molecular dynamics (MDs) simulations, the inhibitory ability of NCI Diversity Set VII natural compounds to VP39 protein was investigated. It should be noted that the computed binding free energy of ligand via molecular docking and linear interaction energy (LIE) approaches are in good agreement with the corresponding experiments with coefficients of R = 0.72 and 0.75, respectively. NSC 319990, NSC 196515 and NSC 376254 compounds were demonstrated that can inhibit MPVX methyltransferase VP39 protein with the similar affinity compared to available inhibitor sinefungin. Moreover, nine residues involving Gln39, Gly68, Gly72, Asp95, Arg97, Val116, Asp138, Arg140 and Asn156 may be argued that they play an important role in binding process of inhibitors to VP39. Communicated by Ramaswamy H. Sarma

Network pharmacology, molecular simulation, and binding free energy calculation-based investigation of Neosetophomone B revealed key targets for the treatment of cancer.

In the current study, Neosetophomone B (NSP-B) was investigated for its anti-cancerous potential using network pharmacology, quantum polarized ligand docking, molecular simulation, and binding free energy calculation. Using SwissTarget prediction, and Superpred, the molecular targets for NSP-B were predicted while cancer-associated genes were obtained from DisGeNet. Among the total predicted proteins, only 25 were reported to overlap with the disease-associated genes. A protein-protein interaction network was constructed by using Cytoscape and STRING databases. MCODE was used to detect the densely connected subnetworks which revealed three sub-clusters. Cytohubba predicted four targets, i.e., fibroblast growth factor , FGF20, FGF22, and FGF23 as hub genes. Molecular docking of NSP-B based on a quantum-polarized docking approach with FGF6, FGF20, FGF22, and FGF23 revealed stronger interactions with the key hotspot residues. Moreover, molecular simulation revealed a stable dynamic behavior, good structural packing, and residues' flexibility of each complex. Hydrogen bonding in each complex was also observed to be above the minimum. In addition, the binding free energy was calculated using the MM/GBSA (Molecular Mechanics/Generalized Born Surface Area) and MM/PBSA (Molecular Mechanics/Poisson-Boltzmann Surface Area) approaches. The total binding free energy calculated using the MM/GBSA approach revealed values of -36.85kcal/mol for the FGF6-NSP-B complex, -43.87kcal/mol for the FGF20-NSP-B complex, and -37.42kcal/mol for the FGF22-NSP-B complex, and -41.91kcal/mol for the FGF23-NSP-B complex. The total binding free energy calculated using the MM/PBSA approach showed values of -30.05kcal/mol for the FGF6-NSP-B complex, -39.62kcal/mol for the FGF20-NSP-B complex, -34.89kcal/mol for the FGF22-NSP-B complex, and -37.18kcal/mol for the FGF23-NSP-B complex. These findings underscore the promising potential of NSP-B against FGF6, FGF20, FGF22, and FGF23, which are reported to be essential for cancer signaling. These results significantly bolster the potential of NSP-B as a promising candidate for cancer therapy.

Comprehensive, Open-Source, and Automated Workflow for Multisite λ-Dynamics in Lead Optimization.

Multisite λ-dynamics (MSLD) is a highly efficient binding free energy calculation method that samples multiple ligands in a single round by assigning different λ values to the alchemical part of each ligand. This method holds great promise for lead optimization (LO) in drug discovery. However, the complex data preparation and simulation process limits its widespread application in diverse protein-ligand systems. To address this challenge, we developed a comprehensive, open-source, and automated workflow for MSLD calculations based on the BLaDE dynamics engine. This workflow incorporates the Ligand Internal and Cartesian coordinate reconstruction-based alignment algorithm (LIC-align) and an optimized maximum common substructure (MCS) search algorithm to accurately generate MSLD multiple topologies with ideal perturbation patterns. Furthermore, our workflow is highly modularized, allowing straightforward integration and extension of various simulation techniques, and is highly accessible to nonexperts. This workflow was validated by calculating the relative binding free energies of large-scale congeneric ligands, many of which have large perturbing groups. The agreement between the calculations and experiments was excellent, with an average unsigned error of 1.08 ± 0.47 kcal/mol. More than 57.1% of the ligands had an error of less than 1.0 kcal/mol, and the perturbations of 6 targets were fully connected via the calculations, while those of 2 targets were connected via both calculations and experimental data. The Pearson correlation coefficient reached 0.88, indicating that the MSLD workflow provides accurate predictions that can guide lead optimization in drug discovery. We also examined the impact of single-site versus multisite perturbations, ligand grouping by perturbing group size, and the position of the anchor atom on the MSLD performance. By integrating our proposed LIC-align and optimized MCS search algorithm along with the coping strategies to handle challenging molecular substructures, our workflow can handle many realistic scenarios more reasonably than all previously published methods. Moreover, we observed that our MSLD workflow achieved similar accuracy to free energy perturbation (FEP) while improving computational efficiency by over 1 order of magnitude in speedup. These findings provide valuable insights and strategies for further MSLD development, making MSLD a competitive tool for lead optimization.

A High-Quality Data Set of Protein-Ligand Binding Interactions Via Comparative Complex Structure Modeling.

High-quality protein-ligand complex structures provide the basis for understanding the nature of noncovalent binding interactions at the atomic level and enable structure-based drug design. However, experimentally determined complex structures are scarce compared with the vast chemical space. In this study, we addressed this issue by constructing the BindingNet data set via comparative complex structure modeling, which contains 69,816 modeled high-quality protein-ligand complex structures with experimental binding affinity data. BindingNet provides valuable insights into investigating protein-ligand interactions, allowing visual inspection and interpretation of structural analogues' structure-activity relationships. It can also be used for evaluating machine-learning-based scoring functions. Our results indicate that machine learning models trained on BindingNet could reduce the bias caused by buried solvent-accessible surface area, as we previously found for models trained on the PDBbind data set. We also discussed strategies to improve BindingNet and its potential utilization for benchmarking the molecular docking methods and ligand binding free energy calculation approaches. The BindingNet complements PDBbind in constructing a sufficient and unbiased protein-ligand binding data set and is freely available at http://bindingnet.huanglab.org.cn.

Control of polymer-protein interactions by tuning the composition and length of polymer chains.

Nanomoduling the 3D shape and chemical functionalities in a synthetic polymer may create recognition cavities for biomacromolecule binding, which serves as an attractive alternative to natural antibodies with much less cost. To obtain fundamental understanding and predict molecular design rules of the polymer antibody, we analyze the complex structure between the biomarker protein epithelial cell adhesion molecule (EpCAM) and a series of polymer ligands via molecular dynamics (MD) simulations. For monomeric ligands, strong enrichment of aromatic residues in protein binding sites is revealed, in line with the reported observations for natural antibodies. Yet, for linear polymers with a growing degree of polymerization, for the first time, a drastic change is revealed on the type of enriched protein residues and the location of protein binding sites, driven by the increasing steric hindrance effect that makes the adsorption of the polymer in the protein exterior feasible. Varying the polymer length and monomeric composition also significantly affects the ligand binding affinity. Here, we have captured three distinct dependences of the ligand binding free energy on the degree of polymerization: for NIPAm based hydrophilic polymers, TBAm dominated hydrophobic polymers and AAc dominated charged polymers. These results can be rationalized by the complex structure and the composition of protein residues at the binding interface. The entire analysis demonstrates unique binding features for polymer ligands and the possibility to modulate their binding sites and affinity by engineering the polymer structure.

Modulation of the conformational landscape of the PDZ3 domain by perturbation on a distal non-canonical α3 helix: decoding the microscopic mechanism of allostery in the PDZ3 domain.

While allosteric signal transduction is crucial for protein signaling and regulation, the dynamic process of allosteric communication remains poorly understood. The third PDZ domain (PDZ stands for the common structural domain shared by the postsynaptic density protein (PSD95), Drosophila disc large tumor suppressor (DlgA), and zonula occludens-1 protein (ZO-1)) serves as a classic example of a single-domain allosteric protein, demonstrating a long-range coupling between the C-terminal α helix (known as the α3 helix) and ligand binding. A molecular level understanding of how the α3 helix modulates the ligand binding affinity of the PDZ3 domain is still lacking. In this study, extensive molecular dynamics simulations corroborated with principal component analysis (PCA), ligand binding free energy calculations, energetic frustration analysis and Markov state model analysis are employed to uncover such molecular details. We demonstrate the definite presence of a binding competent closed-like state in the conformational landscape of wild-type PDZ3. The population modulations of this closed state and other binding incompetent states in the landscape due to α3-truncation/mutation of PDZ3 are explored. A correlation between the closed state population and calculated binding free energy is established, which supports the conformation selection mechanism. Covariance analysis identified the presence of correlated motion between two distant loops (β1-β2 and β2-β3) in the wild-type PDZ3 system, which weakened due to truncation/mutation in the distant α3 helix. It has also been observed that whenever the α3 helix was perturbed, the β2-β3 loop got further away from the binding groove and it is found to be correlated with the binding free energy values. Energetic frustration analysis of the PDZ3 domain also showed that the β2-β3 loop is highly frustrated. Finally, MSM analysis revealed a relevant timescale (closed to open state transition), which is similar to the observed experimental signal transduction timescale for the system. These observations led to the conclusion that the distantly located α3 helix plays a pivotal role in regulating the conformational landscape of the PDZ3 domain, determining the ligand binding affinity and resulting in allosteric behavior of the domain.

Alchemical approach performance in calculating the ligand-binding free energy.

Alchemical binding free energy calculations are one of the most accurate methods for estimating ligand-binding affinity. Assessing the accuracy of the approach over protein targets is one of the most interesting issues. The free energy difference of binding between a protein and a ligand was calculated via the alchemical approach. The alchemical approach exhibits satisfactory accuracy over four targets, including AmpC beta-lactamase (AmpC); glutamate receptor, ionotropic kainate 1 (GluK1); heat shock protein 90 (Hsp90); and severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) main protease (Mpro). In particular, the correlation coefficients between calculated binding free energies and the respective experiments over four targets range from 0.56 to 0.86. The affinity computed via free energy perturbation (FEP) simulations is overestimated over the experimental value. Particularly, the electrostatic interaction free energy rules the binding process of ligands to AmpC and GluK1. However, the van der Waals (vdW) interaction free energy plays an important role in the ligand-binding processes of HSP90 and SARS-CoV-2 Mpro. The obtained results associate with the hydrophilic or hydrophobic properties of the ligands. This observation may enhance computer-aided drug design.

Introducing the Automated Ligand Searcher.

The Automated Ligand Searcher (ALISE) is designed as an automated computational drug discovery tool. To approximate the binding free energy of ligands to a receptor, ALISE includes a three-stage workflow, with each stage involving an increasingly sophisticated computational method: molecular docking, molecular dynamics, and free energy perturbation, respectively. To narrow the number of potential ligands, poorly performing ligands are gradually segregated out. The performance and usability of ALISE are benchmarked for a case study containing known active ligands and decoys for the HIV protease. The example illustrates that ALISE filters the decoys successfully and demonstrates that the automation, comprehensiveness, and user-friendliness of the software make it a valuable tool for improved and faster drug development workflows.

Combining Force Fields and Neural Networks for an Accurate Representation of Chemically Diverse Molecular Interactions.

A key goal of molecular modeling is the accurate reproduction of the true quantum mechanical potential energy of arbitrary molecular ensembles with a tractable classical approximation. The challenges are that analytical expressions found in general purpose force fields struggle to faithfully represent the intermolecular quantum potential energy surface at close distances and in strong interaction regimes; that the more accurate neural network approximations do not capture crucial physics concepts, e.g., nonadditive inductive contributions and application of electric fields; and that the ultra-accurate narrowly targeted models have difficulty generalizing to the entire chemical space. We therefore designed a hybrid wide-coverage intermolecular interaction model consisting of an analytically polarizable force field combined with a short-range neural network correction for the total intermolecular interaction energy. Here, we describe the methodology and apply the model to accurately determine the properties of water, the free energy of solvation of neutral and charged molecules, and the binding free energy of ligands to proteins. The correction is subtyped for distinct chemical species to match the underlying force field, to segment and reduce the amount of quantum training data, and to increase accuracy and computational speed. For the systems considered, the hybrid ab initio parametrized Hamiltonian reproduces the two-body dimer quantum mechanics (QM) energies to within 0.03 kcal/mol and the nonadditive many-molecule contributions to within 2%. Simulations of molecular systems using this interaction model run at speeds of several nanoseconds per day.

Unlocking the Antibiofilm Potential of Natural Compounds by Targeting the NADH:quinone Oxidoreductase WrbA.

Biofilm-dwelling cells endure adverse conditions, including oxidative imbalances. The NADH:quinone oxidoreductase enzyme WrbA has a crucial role in the mechanism of action of antibiofilm molecules such as ellagic and salicylic acids. This study aimed to exploit the potential of the WrbA scaffold as a valuable target for identifying antibiofilm compounds at non-lethal concentrations. A three-dimensional computational model, based on the published WrbA structure, was used to screen natural compounds from a virtual library of 800,000 compounds. Fisetin, morin, purpurogallin, NZ028, and NZ034, along with the reference compound ellagic acid, were selected. The antibiofilm effect of the molecules was tested at non-lethal concentrations evaluating the cell-adhesion of wild-type and WrbA-deprived Escherichia coli strains through fluorochrome-based microplate assays. It was shown that, except for NZ028, all of the selected molecules exhibited notable antibiofilm effects. Purpurogallin and NZ034 showed excellent antibiofilm performances at the lowest concentration of 0.5 μM, in line with ellagic acid. The observed loss of activity and the level of reactive oxygen species in the mutant strain, along with the correlation with terms contributing to the ligand-binding free energy on WrbA, strongly indicates the WrbA-dependency of purpurogallin and NZ034. Overall, the molecular target WrbA was successfully employed to identify active compounds at non-lethal concentrations, thus revealing, for the first time, the antibiofilm efficacy of purpurogallin and NZ034.

Bioassay-guided isolation and in silico study of antihyperlipidemic compounds from Onosma hispidum Wall

Isolation of bioactive compounds from plants and their therapeutic evaluation is crucial in the pursuit of novel phytochemicals and contributes an indispensable role in drug discovery and design. The literature has documented the hypolipidemic effect of numerous Onosma species. Taking that into consideration, the current study was designed to isolate, purify and evaluate the antihyperlipidemic potential of leaves of Onosma hispidum Wall. For the first time, the bioassay-guided isolation led to the separation of 3 compounds that were identified by spectroscopic techniques as o-phthalic acid bis-(2-ethyl decyl)-ester (1), bis (2-ethyloctyl) phthalate (2), and 1,2 benzenedicarboxylic acid bis(2-methyl heptyl) ester (3). Lipase inhibition assay was utilized to scrutinize the antihyperlipidemic potential of methanolic extract fractions and subsequently isolated compounds. Further, the isolated compounds were employed for in silico studies via molecular docking, molecular mechanics with generalized born and surface area solvation (MM-GBSA), and MD simulations with Pancreatic Lipase Colipase (PDB ID: 1LPB). Molecular docking and MM-GBSA of isolated compounds were employed to explain the mode of binding between the protein–ligand complex and binding free energy calculation, respectively. Since compound (3) displayed the best docking score of −6.689 kcal/mol as compared to orlistat −5.529 kcal/mol with PDB: 1LPB. So, it was chosen for MD simulations to evaluate ligand stability and flexibility of the complex which was validated by the fluctuation of α-carbon chain, root mean square deviation (RMSD), root mean square fluctuation (RMSF), and type of interactions involved which authenticated the in vitro lipase inhibitory potential. Overall, in silico and in vitro results validated that compound (3) could be exploited as a promising pancreatic lipase inhibitor. Communicated by Ramaswamy H. Sarma

Enhanced Configurational Sampling Approaches to Alchemical Ligand Binding Free Energy Simulations: Current Status and Challenges.

Ligand binding free energy simulations (LB-FES) have been routine tasks in modern drug discovery campaign. A long-standing challenge for LB-FES is the difficulty in adequately sampling nontrivial environmental reorganizations in response to ligand binding. Therefore, various enhanced configurational sampling (ECS) approaches were devised to speed up fluctuations of relevant slow degrees of freedom (SDOF) and ensure simulation convergence. However, in contrast to the achievements in parametrization, software performance, and workflow automation, efficient ECS methodology suitable for high throughput screening remains in an early stage of development. Here, a review of ECS developments with LB-FES is presented, revisiting current approaches and underlining the major technical pitfalls and challenges. This Perspective focuses on alchemical LB-FES on account of their predominant role in high throughput drug screening as well as the established partnership with ECS. The critical aspects of designing ECS approaches, from both theoretical and applied perspectives, are described. This work is intended to provide a contemporary review of the scientific, technical, and practical issues associated with the accelerating convergence of alchemical LB-FES.

Application and evaluation of molecular docking for aptamer and small molecular interaction - A case study with tetracycline antibiotics

Molecular docking (MD) analysis is currently the most commonly used theoretical simulation method to investigate the interaction of aptamers (receptors) and small molecules (ligands) and understand the recognition mechanism between them at a molecular level. Using the specific aptamers of tetracycline antibiotics (tetracycline (TET), oxytetracycline (OTC), doxycycline (DOC)) as the docking models, three steady-state aptamers of tertiary structures (SATS) were established for each aptamer with the UNAFold and RNAComposer tools. The binding free energy (BFE), docking score (DS), and binding site (base) of the specific ligands (TET, OTC, and DOC) with their respective SATS were obtained by molecular docking. The results revealed one or more binding sites in the established SATS of the aptamers. The BFE and DS of different binding sites of one specific SATS varied significantly. The results also revealed that the site with the highest BFE represented the most dominant binding site, even if it was not the SATS with minimum energy. The BFE values could also be used to evaluate the affinity and specificity of the aptamer to its target. For the first time, this study proposes a method for MD analysis of the aptamer and its target based on different SATS, clarification of the binding mode, and prediction of the binding sites (bases). This study provides a theoretical basis for tailoring; structural optimization; and base modification of aptamers; identifying aptamers with high affinity and specificity.

Long Time Scale Ensemble Methods in Molecular Dynamics: Ligand-Protein Interactions and Allostery in SARS-CoV-2 Targets.

We subject a series of five protein-ligand systems which contain important SARS-CoV-2 targets, 3-chymotrypsin-like protease (3CLPro), papain-like protease, and adenosine ribose phosphatase, to long time scale and adaptive sampling molecular dynamics simulations. By performing ensembles of ten or twelve 10 μs simulations for each system, we accurately and reproducibly determine ligand binding sites, both crystallographically resolved and otherwise, thereby discovering binding sites that can be exploited for drug discovery. We also report robust, ensemble-based observation of conformational changes that occur at the main binding site of 3CLPro due to the presence of another ligand at an allosteric binding site explaining the underlying cascade of events responsible for its inhibitory effect. Using our simulations, we have discovered a novel allosteric mechanism of inhibition for a ligand known to bind only at the substrate binding site. Due to the chaotic nature of molecular dynamics trajectories, regardless of their temporal duration individual trajectories do not allow for accurate or reproducible elucidation of macroscopic expectation values. Unprecedentedly at this time scale, we compare the statistical distribution of protein-ligand contact frequencies for these ten/twelve 10 μs trajectories and find that over 90% of trajectories have significantly different contact frequency distributions. Furthermore, using a direct binding free energy calculation protocol, we determine the ligand binding free energies for each of the identified sites using long time scale simulations. The free energies differ by 0.77 to 7.26 kcal/mol across individual trajectories depending on the binding site and the system. We show that, although this is the standard way such quantities are currently reported at long time scale, individual simulations do not yield reliable free energies. Ensembles of independent trajectories are necessary to overcome the aleatoric uncertainty in order to obtain statistically meaningful and reproducible results. Finally, we compare the application of different free energy methods to these systems and discuss their advantages and disadvantages. Our findings here are generally applicable to all molecular dynamics based applications and not confined to the free energy methods used in this study.

To Design Scalable Free Energy Perturbation Networks, Optimal Is Not Enough.

Drug discovery is accelerated with computational methods such as alchemical simulations to estimate ligand affinities. In particular, relative binding free energy (RBFE) simulations are beneficial for lead optimization. To use RBFE simulations to compare prospective ligands in silico, researchers first plan the simulation experiment, using graphs where nodes represent ligands and graph edges represent alchemical transformations between ligands. Recent work demonstrated that optimizing the statistical architecture of these perturbation graphs improves the accuracy of the predicted changes in the free energy of ligand binding. Therefore, to improve the success rate of computational drug discovery, we present the open-source software package High Information Mapper (HiMap)─a new take on its predecessor, Lead Optimization Mapper (LOMAP). HiMap removes heuristics decisions from design selection and instead finds statistically optimal graphs over ligands clustered with machine learning. Beyond optimal design generation, we present theoretical insights for designing alchemical perturbation maps. Some of these results include that for n number of nodes, the precision of perturbation maps is stable at n·ln(n) edges. This result indicates that even an "optimal" graph can result in unexpectedly high errors if a plan includes too few alchemical transformations for the given number of ligands and edges. And, as a study compares more ligands, the performance of even optimal graphs will deteriorate with linear scaling of the edge count. In this sense, ensuring an A- or D-optimal topology is not enough to produce robust errors. We additionally find that optimal designs will converge more rapidly than radial and LOMAP designs. Moreover, we derive bounds for how clustering reduces cost for designs with a constant expected relative error per cluster, invariant of the size of the design. These results inform how to best design perturbation maps for computational drug discovery and have broader implications for experimental design.

Ligand binding free energy evaluation by Monte Carlo Recursion

The correct evaluation of ligand binding free energies by computational methods is still a very challenging active area of research. The most employed methods for these calculations can be roughly classified into four groups: (i) the fastest and less accurate methods, such as molecular docking, designed to sample a large number of molecules and rapidly rank them according to the potential binding energy; (ii) the second class of methods use a thermodynamic ensemble, typically generated by molecular dynamics, to analyze the endpoints of the thermodynamic cycle for binding and extract differences, in the so-called ‘end-point’ methods; (iii) the third class of methods is based on the Zwanzig relationship and computes the free energy difference after a chemical change of the system (alchemical methods); and (iv) methods based on biased simulations, such as metadynamics, for example. These methods require increased computational power and as expected, result in increased accuracy for the determination of the strength of binding. Here, we describe an intermediate approach, based on the Monte Carlo Recursion (MCR) method first developed by Harold Scheraga. In this method, the system is sampled at increasing effective temperatures, and the free energy of the system is assessed from a series of terms W(b,T), computed from Monte Carlo (MC) averages at each iteration. We show the application of the MCR for ligand binding with datasets of guest-hosts systems (N = 75) and we observed that a good correlation is obtained between experimental data and the binding energies computed with MCR. We also compared the experimental data with an end-point calculation from equilibrium Monte Carlo calculations that allowed us to conclude that the lower-energy (lower-temperature) terms in the calculation are the most relevant to the estimation of the binding energies, resulting in similar correlations between MCR and MC data and the experimental values. On the other hand, the MCR method provides a reasonable view of the binding energy funnel, with possible connections with the ligand binding kinetics, as well. The codes developed for this analysis are publicly available on GitHub as a part of the LiBELa/MCLiBELa project (https://github.com/alessandronascimento/LiBELa).