Binding Energy Distribution Analysis Method: Hamiltonian Replica Exchange with Torsional Flattening for Binding Mode Prediction and Binding Free Energy Estimation

Computational Study on the Wealth Conformations of Hepatitis C Virus Envelope Glycoprotein 2 and Role of Its Glycan Coat

Hepatitis C virus (HCV) is a major cause of chronic liver disease and hepatocellular carcinoma. Antibody development efforts mainly revolve around HCV envelope glycoprotein 2 (E2), which mediates host cell entry by interacting with several cell surface receptors, including CD81. We still have limited knowledge about the structural ensembles and the dynamic behavior of both the CD81 binding sites and the glycans on E2. Here, multiple microsecond-long, all-atom molecular dynamics (MD) simulations, as well as a Markov state model (MSM), were performed to provide an atomistic perspective on the dynamic nature of E2 and its glycans. End-to-end accessibility analyses outline a complete overview of the vulnerabilities of the glycan shield of E2, which may be exploited in therapeutic efforts. Additionally, the Markov state model built from the simulation maps four metastable states for AS412 and three metastable states for the front layer in CD81 binding sites, while binding with HEPC3 would induce a conformation selection for both of them. Overall, this work presents hitherto unseen functional and structural insights into E2 and its glycan coat, providing a new theoretical foundation to control the conformational plasticity of E2 that could be harnessed for vaccine development.

Alchemical Enhanced Sampling with Optimized Phase Space Overlap

An alchemical enhanced sampling (ACES) method has recently been introduced to facilitate importance sampling in free energy simulations. The method achieves enhanced sampling from Hamiltonian replica exchange within a dual topology framework while utilizing new smoothstep softcore potentials. A common sampling problem encountered in lead optimization is the functionalization of aromatic rings that exhibit distinct conformational preferences when interacting with the protein. It is difficult to converge the distribution of ring conformations due to the long time scale of ring flipping events; however, the ACES method addresses this issue by modeling the syn and anti ring conformations within a dual topology. ACES thereby samples the conformer distributions by alchemically tunneling between states, as opposed to traversing a physical pathway with a high rotational barrier. We demonstrate the use of ACES to overcome conformational sampling issues involving ring flipping in ML300-derived noncovalent inhibitors of SARS-CoV-2 Main Protease (Mpro). The demonstrations explore how the use of replica exchange and the choice of softcore selection affects the convergence of the ring conformation distributions. Furthermore, we examine how the accuracy of the calculated free energies is affected by the degree of phase space overlap (PSO) between adjacent states (i.e., between neighboring λ-windows) and the Hamiltonian replica exchange acceptance ratios. Both of these factors are sensitive to the spacing between the intermediate states. We introduce a new method for choosing a schedule of λ values. The method analyzes short "burn-in" simulations to construct a 2D map of the nonlocal PSO. The schedule is obtained by optimizing an alchemical pathway on the 2D map that equalizes the PSO between the λ intervals. The optimized phase space overlap λ-spacing method (Opt-PSO) leads to more numerous end-to-end single passes and round trips due to the correlation between PSO and Hamiltonian replica exchange acceptance ratios. The improved exchange statistics enhance the efficiency of ACES method. The method has been implemented into the FE-ToolKit software package, which is freely available.

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.

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.

Structure-based virtual screening workflow to identify antivirals targeting HIV-1 capsid

We have identified novel HIV-1 capsid inhibitors targeting the PF74 binding site. Acting as the building block of the HIV-1 capsid core, the HIV-1 capsid protein plays an important role in the viral life cycle and is an attractive target for antiviral development. A structure-based virtual screening workflow for hit identification was employed, which includes docking 1.6 million commercially-available drug-like compounds from the ZINC database to the capsid dimer, followed by applying two absolute binding free energy (ABFE) filters on the 500 top-ranked molecules from docking. The first employs the Binding Energy Distribution Analysis Method (BEDAM) in implicit solvent. The top-ranked compounds are then refined using the Double Decoupling method in explicit solvent. Both docking and BEDAM refinement were carried out on the IBM World Community Grid as part of the FightAIDS@Home project. Using this virtual screening workflow, we identified 24 molecules with calculated binding free energies between − 6 and − 12 kcal/mol. We performed thermal shift assays on these molecules to examine their potential effects on the stability of HIV-1 capsid hexamer and found that two compounds, ZINC520357473 and ZINC4119064 increased the melting point of the latter by 14.8 °C and 33 °C, respectively. These results support the conclusion that the two ZINC compounds are primary hits targeting the capsid dimer interface. Our simulations also suggest that the two hit molecules may bind at the capsid dimer interface by occupying a new sub-pocket that has not been exploited by existing CA inhibitors. The possible causes for why other top-scored compounds suggested by ABFE filters failed to show measurable activity are discussed.

Massive-Scale Binding Free Energy Simulations of HIV Integrase Complexes Using Asynchronous Replica Exchange Framework Implemented on the IBM WCG Distributed Network

To perform massive-scale replica exchange molecular dynamics (REMD) simulations for calculating binding free energies of protein-ligand complexes, we implemented the asynchronous replica exchange (AsyncRE) framework of the binding energy distribution analysis method (BEDAM) in implicit solvent on the IBM World Community Grid (WCG) and optimized the simulation parameters to reduce the overhead and improve the prediction power of the WCG AsyncRE simulations. We also performed the first massive-scale binding free energy calculations using the WCG distributed computing grid and 301 ligands from the SAMPL4 challenge for large-scale binding free energy predictions of HIV-1 integrase complexes. In total there are ∼10000 simulated complexes, ∼1 million replicas, and ∼2000 μs of aggregated MD simulations. Running AsyncRE MD simulations on the WCG requires accepting a trade-off between the number of replicas that can be run (breadth) and the number of full RE cycles that can be completed per replica (depth). As compared with synchronous Replica Exchange (SyncRE) running on tightly coupled clusters like XSEDE, on the WCG many more replicas can be launched simultaneously on heterogeneous distributed hardware, but each full RE cycle requires more overhead. We compared the WCG results with that from AutoDock and more advanced RE simulations including the use of flattening potentials to accelerate sampling of selected degrees of freedom of ligands and/or receptors related to slow dynamics due to high energy barriers. We propose a suitable strategy of RE simulations to refine high throughput docking results which can be matched to corresponding computing resources: from HPC clusters, to small or medium-size distributed campus grids, and finally to massive-scale computing networks including millions of CPUs like the resources available on the WCG.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

Protein modeling, molecular network and molecular dynamics study of newly sequenced interleukin-18 (IL-18) gene in Mus musculus

Interleukin-18 (IL-18) belongs to the superfamily of IL-1 protein and exerts a pleiotropic pro-inflammatory effect on the body. Generally, this protein is significantly involved in immune defense during infection in cells, but sometimes its anomalous activities produce some inflammatory diseases like rheumatoid arthritis and Crohn's disease. In the present study, the IL-18 gene was isolated from mice and was subsequently cloned and sequenced. Further, the network analysis was carried out to explore the functional role of IL-18 protein in animals. The 3D protein structure of the IL-18 protein was generated and docked with appropriate 3-([3-cholamidopropyl]dimethylammonio)-1-propanesulfonate (CPS) ligand. Later the complex structure of the protein was subjected to molecular dynamics simulation (MDS) for 50 ns to determine the effect of ligand on protein. The network analysis explored the correlation of IL-18 protein with others proteins and their involvement in the different significant pathway to defend the cell from various diseases. As confirmed by MDS, the CPS:IL-18 complex was found to be highly stable. Our results further indicated that CPS ligand has the potential to act as a drug molecule, in future, for counteracting IL-18 activity. To date, no structural details were available for animal IL-18. Hence, the finding of this study will be useful in broadening the horizon towards a better understanding of the functional and structural aspects of IL-18 in animals.

Reduced Free Energy Perturbation/Hamiltonian Replica Exchange Molecular Dynamics Method with Unbiased Alchemical Thermodynamic Axis

Replica-exchange molecular dynamics (REMD) has been proven to efficiently improve the convergence of free-energy perturbation (FEP) calculations involving considerable reorganization of their surrounding. We previously introduced the FEP/(λ,H)-REMD algorithm for ligand binding, in which replicas along the alchemical thermodynamic coupling axis λ were expanded as a series of Hamiltonian boosted replicas along a second axis to form a two-dimensional replica-exchange exchange map [Jiang, W.; Roux, B., J. Chem. Theory Comput. 2010, 6 (9), 2559-2565]. Aiming to achieve a similar performance at a lower computational cost, we propose here a modified version of this algorithm in which only the end-states along the alchemical axis are augmented by boosted replicas. The reduced FEP/(λ,H)-REMD method with one-dimensional unbiased alchemical thermodynamic coupling axis λ is implemented on the basis of generic multiple copy algorithm (MCA) module of the biomolecular simulation program NAMD. The flexible MCA framework of NAMD enables a user to design customized replica-exchange patterns through Tcl scripting in the context of a highly parallelized simulation program without touching the source code. Two Hamiltonian tempering boosting scheme were examined with the new algorithm: a first one based on potential energy rescaling of a preidentified "solute" and a second one via the introduction of flattening torsional free-energy barriers. As two illustrative examples with reliable experiment data, the absolute binding free energies of p-xylene and n-butylbenzene to the nonpolar cavity of the L99A mutant of T4 lysozyme were calculated. The tests demonstrate that the new protocol efficiently enhances the sampling of torsional motions for backbone and side chains around the binding pocket and accelerates the convergence of the free-energy computations.

Improving Prediction Accuracy of Binding Free Energies and Poses of HIV Integrase Complexes Using the Binding Energy Distribution Analysis Method with Flattening Potentials

To accelerate conformation sampling of slow dynamics from receptor or ligand, we introduced flattening potentials on selected bonded and nonbonded intramolecular interactions to the binding energy distribution analysis method (BEDAM) for calculating absolute binding free energies of protein-ligand complexes using an implicit solvent model and implemented flattening BEDAM using the asynchronous replica exchange (AsyncRE) framework for performing large scale replica exchange molecular dynamics (REMD) simulations. The advantage of using the flattening feature to reduce high energy barriers was exhibited first by the p-xylene-T4 lysozyme complex, where the intramolecular interactions of a protein side chain on the binding site were flattened to accelerate the conformational transition of the side chain from the trans to the gauche state when the p-xylene ligand is present in the binding site. Much more extensive flattening BEDAM simulations were performed for 53 experimental binders and 248 nonbinders of HIV-1 integrase which formed the SAMPL4 challenge, with the total simulation time of 24.3 μs. We demonstrated that the flattening BEDAM simulations not only substantially increase the number of true positives (and reduce false negatives) but also improve the prediction accuracy of binding poses of experimental binders. Furthermore, the values of area under the curve (AUC) of receiver operating characteristic (ROC) and the enrichment factors at 20% cutoff calculated from the flattening BEDAM simulations were improved significantly in comparison with that of simulations without flattening as we previously reported for the whole SAMPL4 database. Detailed analysis found that the improved ability to discriminate the binding free energies between the binders and nonbinders is due to the fact that the flattening simulations reduce the reorganization free energy penalties of binders and decrease the overlap of binding free energy distributions of binders relative to that of nonbinders. This happens because the conformational ensemble distributions for both the ligand and protein in solution match those at the fully coupled (complex) state more closely when the systems are more fully sampled after the flattening potentials are applied to the intermediate states.

The effect of β-glucan and its potential analog on the structure of Dectin-1 receptor

Dectin-1 is a recently discovered pattern-recognition receptor that plays an important role in antifungal innate immunity, which acts a specific receptor for β-glucan (BG). The present study, aimed at clarifying effect of BG and a new analog, maltotriose (MT) on Dectin-1 receptor. We implemented molecular docking of MT on Dectin-1 along with model-independent all-atom-molecular dynamics simulations. Simulations were carried out at three levels of complexity: (1) Apo-Dectin-1; (2) BG:Dectin-1; (3) MT:Dectin-1. All three system complexes were undergone stability check before showing a comparative analysis. A characteristic feature, noted for the MT:Dectin-1, is a shifting of loops (loop1 and loop2) orientation towards atoms of MT, a broad interaction suggested a robust and tight binding on comparison with BG:Dectin-1. Free energy estimation corroborated the observation, which furthermore, made a close agreement by revealing contribution of energy components of interacting residues. In addition, cluster analysis of complexes exhibit a smooth continuous transition to a new confirmation, represented by a series of clusters each having a longer lifetime. Principal component analysis revealed a broken pipe at binding site of BG:Dectin-1 during movement of atoms whereas in MT:Dectin-1 exhibited wide band and high amplitude motion of atoms in trajectory, was due to loop orientation toward MT. Observation was further shown by measuring distances and hydrogen binding calculation. Simulations of the BG:Dectin-1 and MT:Dectin-1 complex revealed first time the influence of BG and MT ligands. This study might extend the knowledge of the BG and MT interaction on Dectin-1 and proposed further potential bioassay of MT.

Large scale free energy calculations for blind predictions of protein-ligand binding: the D3R Grand Challenge 2015

We describe binding free energy calculations in the D3R Grand Challenge 2015 for blind prediction of the binding affinities of 180 ligands to Hsp90. The present D3R challenge was built around experimental datasets involving Heat shock protein (Hsp) 90, an ATP-dependent molecular chaperone which is an important anticancer drug target. The Hsp90 ATP binding site is known to be a challenging target for accurate calculations of ligand binding affinities because of the ligand-dependent conformational changes in the binding site, the presence of ordered waters and the broad chemical diversity of ligands that can bind at this site. Our primary focus here is to distinguish binders from nonbinders. Large scale absolute binding free energy calculations that cover over 3000 protein-ligand complexes were performed using the BEDAM method starting from docked structures generated by Glide docking. Although the ligand dataset in this study resembles an intermediate to late stage lead optimization project while the BEDAM method is mainly developed for early stage virtual screening of hit molecules, the BEDAM binding free energy scoring has resulted in a moderate enrichment of ligand screening against this challenging drug target. Results show that, using a statistical mechanics based free energy method like BEDAM starting from docked poses offers better enrichment than classical docking scoring functions and rescoring methods like Prime MM-GBSA for the Hsp90 data set in this blind challenge. Importantly, among the three methods tested here, only the mean value of the BEDAM binding free energy scores is able to separate the large group of binders from the small group of nonbinders with a gap of 2.4 kcal/mol. None of the three methods that we have tested provided accurate ranking of the affinities of the 147 active compounds. We discuss the possible sources of errors in the binding free energy calculations. The study suggests that BEDAM can be used strategically to discriminate binders from nonbinders in virtual screening and to more accurately predict the ligand binding modes prior to the more computationally expensive FEP calculations of binding affinity.