Ligand Binding Free Energy Research Articles

Estimating the ligand-binding affinity via λ-dependent umbrella sampling simulations.

The umbrella sampling (US) approach has been demonstrated to be a very efficient method for estimating the ligand-binding affinity. However, most of the calculated values overestimate experimental ones that are probably caused by the inaccurate representation of the interaction between the ligand and the surrounding molecules. The issue can be resolved via the implementation aspects of λ-alteration simulation into the US approach, which we call the λ-dependent umbrella sampling (λUS) scheme. In particular, the electrostatic and van der Waals interactions were simultaneously changed by using the coupling parameter λ during λUS simulations. The mean value of obtained results, kcal mol-1 , is in good fitting to the mean value of respective experiments, ∆GEXP = - 11.26 ± 0.89 kcal mol-1 . Moreover, the correlation between the proposed approach and experiment is quite good with a value of . The λUS scheme significantly enhances the calculated accuracy since the RMSE of the proposed scheme is smaller than traditional US simulations, kcal mol-1 versus kcal mol-1 . Furthermore, the precision is increased since the computed error via λUS approach, kcal mol-1 , was smaller than those of the US simulation, kcal mol-1 . Overall, the proposed approach perhaps provides an efficient way to accurately and precisely estimate the ligand-binding free energy.

Receptor-Ligand Binding Free Energies from a Consecutive Histograms Monte Carlo Sampling Method.

To obtain accurate and converged free energy calculations for ligand binding to biomolecular systems requires validated force fields and extensive sampling of the energy landscape, which requires exhaustive and effective conformational searching methods. Herein, we introduce the consecutive histograms Monte Carlo (CHMC) sampling protocol that generates receptor-ligand binding modes within a series of continuously distributed sampling units ranging from placement near the geometric center of the receptor's binding site to fully unbound states. This protocol employs independent energy-state sampling for calculating the ensemble energy within every predefined location along the receptor-ligand dissociation pathway, without the need to traverse the energy barriers as in molecular dynamic simulations during the dissociation procedure. We applied this method to a set of selected receptor targets with their corresponding ligands providing detailed studies of molecular binding free energy predictions. The results show that the CHMC gives an excellent accounting of the free energy surfaces and binding free energies at a reasonable computational cost.

Activity of Thiazine substituted 9-anilinoacridines against Corona virus (COVID19): An In-silico approach

Coronavirus Disease 2019 (COVID-19), a life-threatening viral disease affected first in Wuhan, China, and quickly spread to more than 200 countries in the world in the year 2020. So many scientists are trying to discover novel drugs and vaccines for coronavirus and treatment for COVID-19. In the present article, in-silico studies have been performed to explore the binding modes of Thiazine substituted 9-anilinoacridines (1a-z) against SARS CoV 2 main protease (PDB id - 5R82) targeting the coronavirus using Schrodinger suit 2019-4. The molecular docking studies are performed by Glide module, in-silico ADMET screening was performed by Qik prop module, and the binding free energy of ligands was calculated using PRIME MM-GB/SA module of Schrodinger suite 2019-4, Maestro 21.2 version. From the in-silico results, Thiazine substituted 9-anilinoacridines like 1m, 1j, 1s and 1b are significantly active against SARS CoV 2 main protease with Glide score more than -5.4 when compared with the currently recommended drug for COVID19, Hydroxychloroquine (G score -5.47). The docking results of the Thiazine substituted 9-anilinoacridines exhibited similar mode of interactions with COVID19 and the residues GLN19, THR24, THR25, THR26, LEU27, HIE41, SER46, MET49, ASN142, GLN143, HIE164, MET165, ASP187, ARG188 and GLN189, play a crucial role in binding with ligands.

Allosteric Mechanisms of Nonadditive Substituent Contributions to Protein-Ligand Binding

Quantifying chemical substituent contributions to ligand-binding free energies is challenging due to nonadditive effects. Protein allostery is a frequent cause of nonadditivity, but the underlying allosteric mechanisms often remain elusive. Here, we propose a general NMR-based approach to elucidate such mechanisms and we apply it to the HCN4 ion channel, whose cAMP-binding domain is an archetypal conformational switch. Using NMR, we show that nonadditivity arises not only from concerted conformational transitions, but also from conformer-specific effects, such as steric frustration. Our results explain how affinity-reducing functional groups may lead to affinity gains if combined. Surprisingly, our approach also reveals that nonadditivity depends markedly on the receptor conformation. It is negligible for the inhibited state but highly significant for the active state, opening new opportunities to tune potency and agonism of allosteric effectors.

Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD): Characterization of Ligand Binding Thermodynamics and Kinetics.

Calculations of ligand binding free energies and kinetic rates are important for drug design. However, such tasks have proven challenging in computational chemistry and biophysics. To address this challenge, we have developed a new computational method, ligand Gaussian accelerated molecular dynamics (LiGaMD), which selectively boosts the ligand nonbonded interaction potential energy based on the Gaussian accelerated molecular dynamics (GaMD) enhanced sampling technique. Another boost potential could be applied to the remaining potential energy of the entire system in a dual-boost algorithm (LiGaMD_Dual) to facilitate ligand binding. LiGaMD has been demonstrated on host-guest and protein-ligand binding model systems. Repetitive guest binding and unbinding in the β-cyclodextrin host were observed in hundreds-of-nanosecond LiGaMD_Dual simulations. The calculated guest binding free energies agreed excellently with experimental data with <1.0 kcal/mol errors. Compared with converged microsecond-time scale conventional molecular dynamics simulations, the sampling errors of LiGaMD_Dual simulations were also <1.0 kcal/mol. Accelerations of ligand kinetic rate constants in LiGaMD simulations were properly estimated using Kramers' rate theory. Furthermore, LiGaMD allowed us to capture repetitive dissociation and binding of the benzamidine inhibitor in trypsin within 1 μs simulations. The calculated ligand binding free energy and kinetic rate constants compared well with the experimental data. In summary, LiGaMD provides a powerful enhanced sampling approach for characterizing ligand binding thermodynamics and kinetics simultaneously, which is expected to facilitate computer-aided drug design.

Toward Automated Free Energy Calculation with Accelerated Enveloping Distribution Sampling (A-EDS).

Free-energy perturbation (FEP) methods are commonly used in drug design to calculate relative binding free energies of different ligands to a common host protein. Alchemical ligand transformations are usually performed in multiple steps which need to be chosen carefully to ensure sufficient phase-space overlap between neighboring states. With one-step or single-step FEP techniques, a single reference state is designed that samples phase-space not only representative of a full transformation but also ideally resembles multiple ligand end states and hence allows for efficient multistate perturbations. Enveloping distribution sampling (EDS) is one example for such a method in which the reference state is created by a mathematical combination of the different ligand end states based on solid statistical mechanics. We have recently proposed a novel approach to EDS which enables efficient barrier crossing between the different end states, termed accelerated EDS (A-EDS). In this work, we further simplify the parametrization of the A-EDS reference state and demonstrate the automated calculation of multiple free-energy differences between different ligands from a single simulation in three different well-described drug design model systems.

ProtoCaller: Robust Automation of Binding Free Energy Calculations.

ProtoCaller is a Python library distributed through Anaconda which automates relative protein-ligand binding free energy calculations in GROMACS. It links a number of popular specialized tools used to perform protein setup and parametrization, such as PDB2PQR, Modeller, and AmberTools. ProtoCaller supports commonly used AMBER force fields with additional cofactor parameters, and AM1-BCC is used to derive ligand charges. ProtoCaller also comes with an extensive PDB parser, an enhanced maximum common substructure algorithm providing improved ligand-ligand mapping, and a light GROMACS wrapper for running multiple molecular dynamics simulations. ProtoCaller is highly relevant to most researchers in the field of biomolecular simulation, allowing a customizable balance between automation and user intervention.

Thermodynamic contribution of iodine atom to the binding of heterogeneously polyhalogenated benzotriazoles by the catalytic subunit of human protein kinase CK2.

A series of novel benzotriazole derivatives containing iodine atom(s) were synthesized. The binding of these compounds to the catalytic subunit of human protein kinase CK2 was evaluated using differential scanning fluorimetry. The obtained thermodynamic data were compared with those determined previously for the brominated and chlorinated benzotriazole analogues to get a deeper insight into the thermodynamic contribution of iodine substitution to the free energy of ligand binding. We have shown that iodine atom(s) attached to the benzene ring of benzotriazole enhance(s) its binding by the target protein. This effect is the strongest when two iodine atoms are attached at positions peripheral to the triazole ring, which according to the structures deposited in protein data bank may be indicative for the formation of the halogen bond between iodine and carbonyl groups of residues located in the hinge region of the protein. Finally, quantitative structure-activity relationship analysis pointed the solute hydrophobicity as the main factor contributing to the binding affinity.

Elucidating the multiple roles of hydration for accurate protein-ligand binding prediction via deep learning

Accurate and efficient prediction of protein-ligand interactions has been a long-lasting dream of practitioners in drug discovery. The insufficient treatment of hydration is widely recognized to be a major limitation for accurate protein-ligand scoring. Using an integration of molecular dynamics simulations on thousands of protein structures with novel big-data analytics based on convolutional neural networks and deep Taylor decomposition, we consistently identify here three different patterns of hydration to be essential for protein-ligand interactions. In addition to desolvation and water-mediated interactions, the formation of enthalpically favorable networks of first-shell water molecules around solvent-exposed ligand moieties is identified to be essential for protein-ligand binding. Despite being currently neglected in drug discovery, this hydration phenomenon could lead to new avenues in optimizing the free energy of ligand binding. Application of deep neural networks incorporating hydration to docking provides 89% accuracy in binding pose ranking, an essential step for rational structure-based drug design.

Correcting electrostatic artifacts due to net-charge changes in the calculation of ligand binding free energies.

Alchemically derived free energies are artifacted when the perturbed moiety has a nonzero net charge. The source of the artifacts lies in the effective treatment of the electrostatic interactions within and between the perturbed atoms and remaining (partial) charges in the simulated system. To treat the electrostatic interactions effectively, lattice-summation (LS) methods or cutoff schemes in combination with a reaction-field contribution are usually employed. Both methods render the charging component of the calculated free energies sensitive to essential parameters of the system like the cutoff radius or the box side lengths. Here, we discuss the results of three previously published studies of ligand binding. These studies presented estimates of binding free energies that were artifacted due to the charged nature of the ligands. We show that the size of the artifacts can be efficiently calculated and raw simulation data can be corrected. We compare the corrected results with experimental estimates and nonartifacted estimates from path-sampling methods. Although the employed correction scheme involves computationally demanding continuum-electrostatics calculations, we show that the correction estimate can be deduced from a small sample of configurations rather than from the entire ensemble. This observation makes the calculations of correction terms feasible for complex biological systems. To show the general applicability of the proposed procedure, we also present results where the correction scheme was used to correct independent free energies obtained from simulations employing a cutoff scheme or LS electrostatics. In this work, we give practical guidelines on how to apply the appropriate corrections easily.

Ligand binding free energy and kinetics calculation in 2020

AbstractLigand/protein binding (LPB) is a major topic in medicine, chemistry and biology. Since the advent of computers, many scientists have put efforts in developing theoretical models that could decode the alphabet of the LPB interaction. The success of this task passes by the resolution of the molecular mechanism of LPB. In the past century, major attention was dedicated to the thermodynamics of LPB, while more recent studies have revealed that ligand (un)binding kinetics is at least as important as ligand binding thermodynamics in determining the drug in vivo efficacy. In the present review, we introduce the most widely used computational methods to study LPB thermodynamics and kinetics. It is important to say that no method outperforms another, they all have pros and cons and the choice of the user should take carefully into account the system under investigation, the available structural and experimental data, and the goal of the research. A perspective on future directions of method development and research on LPB concludes the discussion.This article is categorized under: Molecular and Statistical Mechanics &gt; Free Energy Methods Structure and Mechanism &gt; Computational Biochemistry and Biophysics Molecular and Statistical Mechanics &gt; Molecular Dynamics and Monte‐Carlo Methods

Taking Water into Account with the Fragment Molecular Orbital Method.

This chapter describes the current status of development of the fragment molecular orbital (FMO) method for analyzing the electronic state and intermolecular interactions of biomolecular systems in solvent. The orbital energies and the inter-fragment interaction energies (IFIEs) for a specific molecular structure can be obtained directly by performing FMO calculations by exposing water molecules and counterions around biomolecular systems. Then, it is necessary to pay attention to the thickness of the water shell surrounding the biomolecules. The single-point calculation for snapshots from MD trajectory does not incorporate the effects of temperature and configurational fluctuation, but the SCIFIE (statistically corrected IFIE) method is proposed as a many-body correlated method that partially compensates for this deficiency. Furthermore, implicit continuous dielectric models have been developed as effective approaches to incorporating the screening effect of the solvent in thermal equilibrium, and we illustrate their usefulness for theoretical evaluation of IFIEs and ligand-binding free energy on the basis of the FMO-PBSA (Poisson-Boltzmann surface area) method and other computational methods.

Modelling flexible protein-ligand binding in p38α MAP kinase using the QUBE force field.

The quantum mechanical bespoke (QUBE) force field is used to retrospectively calculate the relative binding free energies of a series of 17 flexible inhibitors of p38α MAP kinase. The size and flexibility of the chosen molecules represent a stringent test of the derivation of force field parameters from quantum mechanics, and enhanced sampling is required to reduce the dependence of the results on the starting structure. Competitive accuracy with a widely-used biological force field is achieved, indicating that quantum mechanics derived force fields are approaching the accuracy required to provide guidance in prospective drug discovery campaigns.

Oversampling Free Energy Perturbation Simulation in Determination of the Ligand-Binding Free Energy.

Determination of the ligand-binding affinity is an extremely interesting problem. Normally, the free energy perturbation (FEP) method provides an appropriate result. However, it is of great interest to improve the accuracy and precision of this method. In this context, temperature replica exchange molecular dynamics implementation of the FEP computational approach, which we call replica exchange free energy perturbation (REP) was proposed. In particular, during REP simulations, the system can easily escape from being trapped in local minima by exchanging configurations with high temperatures, resulting in significant improvement in the accuracy and precision of protein-ligand binding affinity calculations. The distribution of the decoupling free energy was enlarged, and its mean values were decreased. This results in changes in the magnitude of the calculated binding free energies as well as in alteration in the binding mechanism. Moreover, the REP correlation coefficient with respect to experiment ( RREP = 0.85 ± 0.15) is significantly boosted in comparison with the FEP one ( RFEP = 0.64 ± 0.30). Furthermore, the root-mean-square error (RMSE) of REP is also smaller than FEP, RMSEREP = 4.28 ± 0.69 versus RMSEFEP = 5.80 ± 1.11 kcal/mol, respectively. © 2019 Wiley Periodicals, Inc.

Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity

Molecular docking simulation is a useful tool in the prediction of protein-ligand binding affinity on a large scale and has great potential in various application fields such as virtual screening of potential drug molecules. However, the reliability of molecular docking is still weak in the estimation of ligand-binding free energy, which limits the applicability of molecular docking simulation. Ligand torsion number is related to the flexibility of ligand and generally incorporated as a crucial variable in the thermodynamic function of binding free energy. In this study, we investigated how the ligand torsion number has influence on the binding affinity prediction of AutoDock, a popular molecular docking simulation tool. The pKd values of various protein-ligands were estimated by using the binding free energy function of AutoDock and compared with their experimental pKd values. The torsion number dependent comparison showed that the predicted binding affinities were mostly underestimated in the complexes of higher torsion numbers, whereas the underestimated and overestimated cases were relatively balanced at relatively lower torsion numbers. A new weight factor for torsion-free energy term of binding energy function was determined and introduced to make correction to the underestimation of binding affinity of ligands with high torsion numbers. It is expected that the torsion number dependent deviation pattern of AutoDock and its correction strategy are useful in the large-scale validation of protein-ligand binding affinity.

On Restraints in End-Point Protein-Ligand Binding Free Energy Calculations.

The impact of harmonic restraints on protein heavy atoms and ligand atoms on end-point free energy calculations is systematically characterized for 54 protein-ligand complexes. We observe that stronger restraints reduce the equilibration time and statistical inefficiency, suppress conformational sampling, influence correlation with experiment, and monotonically decrease the estimated loss of entropy upon binding, leading to stronger estimated binding free energies in most systems. A statistical estimator that reweights for the biasing potential and includes data prior to the estimated equilibration time has the highest correlation with experiment. A spring constant of 20 cal mol-1 Å-2 maintains a near-native energy landscape and suppresses artifactual energy minima while minimally limiting thermal fluctuations about the crystal structure. © 2019 Wiley Periodicals, Inc.

Inhibition Mechanisms of Rhodococcus Erythropolis 2'-Hydroxybiphenyl-2-sulfinate Desulfinase (DszB).

Naturally occurring enzymatic pathways enable highly specific, rapid thiophenic sulfur cleavage occurring at ambient temperature and pressure, which may be harnessed for the desulfurization of petroleum-based fuel. One pathway found in bacteria is a four-step catabolic pathway (the 4S pathway) converting dibenzothiophene (DBT), a common crude oil contaminant, into 2-hydroxybiphenyl (HBP) without disrupting the carbon-carbon bonds. 2'-Hydroxybiphenyl-2-sulfinate desulfinase (DszB), the rate-limiting enzyme in the enzyme cascade, is capable of selectively cleaving carbon-sulfur bonds. Accordingly, understanding the molecular mechanisms of DszB activity may enable development of the cascade as industrial biotechnology. Based on crystallographic evidence, we hypothesized that DszB undergoes an active site conformational change associated with the catalytic mechanism. Moreover, we anticipated this conformational change is responsible, in part, for enhancing product inhibition. Rhodococcus erythropolis IGTS8 DszB was recombinantly produced and purified via Escherichia coli BL21 to test these hypotheses. Activity and the resulting conformational change of DszB in the presence of HBP were evaluated. The activity of recombinant DszB was comparable to the natively expressed enzyme and was inhibited via competitive binding of the product, HBP. Using circular dichroism, global changes in DszB conformation were monitored in response to HBP concentration, which indicated that both product and substrate produced similar structural changes. Molecular dynamics (MD) simulations and free energy perturbation with Hamiltonian replica exchange molecular dynamics (FEP/λ-REMD) calculations were used to investigate the molecular-level phenomena underlying the connection between conformation change and kinetic inhibition. In addition to the HBP, MD simulations of DszB bound to common, yet structurally diverse, crude oil contaminants 2',2-biphenol (BIPH), 1,8-naphthosultam (NTAM), 2-biphenyl carboxylic acid (BCA), and 1,8-naphthosultone (NAPO) were performed. Analysis of the simulation trajectories, including root-mean-square fluctuation (RMSF), center of mass (COM) distances, and strength of nonbonded interactions, when compared with FEP/λ-REMD calculations of ligand binding free energy, showed excellent agreement with experimentally determined inhibition constants. Together, the results show that the combination of a molecule's hydrophobicity and nonspecific interactions with nearby functional groups contributes to a competitive inhibition mechanism that locks DszB in a closed conformation and precludes substrate access to the active site.

Development of a Robust Indirect Approach for MM → QM Free Energy Calculations That Combines Force-Matched Reference Potential and Bennett's Acceptance Ratio Methods.

We use the PBE0/6-31G* density functional method to perform ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations under periodic boundary conditions with rigorous electrostatics using the ambient potential composite Ewald method in order to test the convergence of MM → QM/MM free energy corrections for the prediction of 17 small-molecule solvation free energies and eight ligand binding free energies to T4 lysozyme. The "indirect" thermodynamic cycle for calculating free energies is used to explore whether a series of reference potentials improve the statistical quality of the predictions. Specifically, we construct a series of reference potentials that optimize a molecular mechanical (MM) force field's parameters to reproduce the ab initio QM/MM forces from a QM/MM simulation. The optimizations form a systematic progression of successively expanded parameters that include bond, angle, dihedral, and charge parameters. For each reference potential, we calculate benchmark quality reference values for the MM → QM/MM correction by performing the mixed MM and QM/MM Hamiltonians at 11 intermediate states, each for 200 ps. We then compare forward and reverse application of Zwanzig's relation, thermodynamic integration (TI), and Bennett's acceptance ratio (BAR) methods as a function of reference potential, simulation time, and the number of simulated intermediate states. We find that Zwanzig's equation is inadequate unless a large number of intermediate states are explicitly simulated. The TI and BAR mean signed errors are very small even when only the end-state simulations are considered, and the standard deviations of the TI and BAR errors are decreased by choosing a reference potential that optimizes the bond and angle parameters. We find a robust approach for the data sets of fairly rigid molecules considered here is to use bond + angle reference potential together with the end-state-only BAR analysis. This requires QM/MM simulations to be performed in order to generate reference data to parametrize the bond + angle reference potential, and then this same simulation serves a dual purpose as the full QM/MM end state. The convergence of the results with respect to time suggests that computational resources may be used more efficiently by running multiple simulations for no more than 50 ps, rather than running one long simulation.

Does the partial molar volume of a solute reflect the free energy of hydrophobic solvation?

Halogenated heterocyclic ligands are widely used as the potent and frequently selective inhibitors of protein kinases. However, the exact contribution of the hydrophobic solvation of a free ligand is rarely accounted for the balance of interactions contributing to the free energy of ligand binding. Herein, we propose a new experimental method based on volumetric data to estimate the hydrophobicity of a ligand. We have tested this approach for a series of ten variously halogenated benzotriazoles, the binding affinity of which to the target protein kinase CK2 was assessed with the use of thermal shift assay. According to the hierarchical clustering procedure, the excess volume, defined as the difference between the experimentally determined partial molar volume and the calculated in silico molecular volume, was found to be distant from any commonly used hydrophobicity descriptors of the ligand. The excess volume, however, properly predicts solute binding affinity. On the way, we have proved that the binding of halogenated benzotriazoles to the protein kinase CK2 is driven mostly by hydrophobic interactions.

Tackling Halogenated Species with PBSA: Effect of Emulating the σ-Hole.

To model halogen-bond phenomena using classical force fields, an extra point (EP) of charge is frequently introduced at a given distance from the halogen (X) to emulate the σ-hole. The resulting molecular dynamics (MD) trajectories can be used in subsequent molecular mechanics (MM) combined with Poisson-Boltzmann and surface area calculations (PBSA) to estimate protein-ligand binding free energies (Δ Gbind). While EP addition improves the MM/MD description of halogen-containing systems, its effect on the calculation of solvation free energies (Δ Gsolv) using the PBSA approach is yet to be assessed. As the PBSA calculations depend, among other parameters, on the empirical assignment of radii (PB radii), a problematic issue arises, since standard halogen radii are smaller than the typical X···EP distances, thus placing the EP within the solvent dielectric. Herein, we took a common literature EP parametrization scheme, which uses X···EP = Rmin and RESP charges in the context of GAFF, and performed a comprehensive study on the performance of PBSA (using three different setups) in the calculation of Δ Gsolv values for 142 halogenated compounds (bearing Cl, Br, or I) for which the experimental values are known. By conducting an optimization (minimizing the error against experimental values), we provide a new optimized set of halogen PB radii, for each PBSA setup, that should be used in the context of the aforementioned scenario. A simultaneous optimization of PB radii and X···EP distances shows that a wide range of distance/radius pairs can be used without significant loss of accuracy, therefore laying the basis for expanding this halogen radii optimization strategy to other force fields and EP implementations. As ligand Δ Gsolv estimation is an important term in the determination of protein-ligand Δ Gbind, this work is particularly relevant in the framework of structure-based virtual screening and related computer-aided drug design routines.