Fast Equilibration of Water between Buried Sites and the Bulk by Molecular Dynamics with Parallel Monte Carlo Water Moves on Graphical Processing Units

Impact of Varying Velocities and Solvation Boxes on Alchemical Free-Energy Simulations

Alchemical free-energy perturbation (FEP) is an accurate and thermodynamically stringent way to estimate relative energies for the binding of small ligands to biological macromolecules. It has repeatedly been pointed out that a single simulation normally stays near the starting point in phase space and therefore underestimates the uncertainty of the results. Therefore, it is better to run an ensemble of independent simulations. Traditionally, such an ensemble has been generated by using different starting velocities. We argue that it is better to use also other random choices made during the setup of the simulations, in particular the solvation of the solute. We show here that such solvent-induced independent simulations (SIS) sometimes give a larger standard deviation and slightly different results for the binding of 42 ligands to five different proteins, viz. human N-terminal bromodomain 4, the Leu99Ala mutant of T4 lysozyme, dihydrofolate reductase, blood-clotting factor Xa, and ferritin. SIS does not involve any increase in the time consumption. Therefore, we strongly recommend the use of SIS (in addition to different velocities) to start independent simulations. Other random or uncertain choices in the setup of the simulated systems, e.g., the selection of residues with alternative conformations or positions of added protons, may also be used to enhance the variation in independent simulations.

Efficient Sampling of Cavity Hydration in Proteins with Nonequilibrium Grand Canonical Monte Carlo and Polarizable Force Fields

Prediction of the hydration levels of protein cavities and active sites is important to both mechanistic analysis and ligand design. Due to the unique microscopic environment of these buried water molecules, a polarizable model is expected to be crucial for an accurate treatment of protein internal hydration in simulations. Here we adapt a nonequilibrium candidate Monte Carlo approach for conducting grand canonical Monte Carlo simulations with the Drude polarizable force field. The GPU implementation enables the efficient sampling of internal cavity hydration levels in biomolecular systems. We also develop an enhanced sampling approach referred to as B-walking, which satisfies detailed balance and readily combines with grand canonical integration to efficiently calculate quantitative binding free energies of water to protein cavities. Applications of these developments are illustrated in a solvent box and the polar ligand binding site in trypsin. Our simulation results show that including electronic polarization leads to a modest but clear improvement in the description of water position and occupancy compared to the crystal structure. The B-walking approach enhances the range of water sampling in different chemical potential windows and thus improves the accuracy of water binding free energy calculations.

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.

Accelerated Enveloping Distribution Sampling to Probe the Presence of Water Molecules

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

WaterKit: Thermodynamic Profiling of Protein Hydration Sites

Water desolvation is one of the key components of the free energy of binding of small molecules to their receptors. Thus, understanding the energetic balance of solvation and desolvation resulting from individual water molecules can be crucial when estimating ligand binding, especially when evaluating different molecules and poses as done in High-Throughput Virtual Screening (HTVS). Over the most recent decades, several methods were developed to tackle this problem, ranging from fast approximate methods (usually empirical functions using either discrete atom-atom pairwise interactions or continuum solvent models) to more computationally expensive and accurate ones, mostly based on Molecular Dynamics (MD) simulations, such as Grid Inhomogeneous Solvation Theory (GIST) or Double Decoupling. On one hand, MD-based methods are prohibitive to use in HTVS to estimate the role of waters on the fly for each ligand. On the other hand, fast and approximate methods show an unsatisfactory level of accuracy, with low agreement with results obtained with the more expensive methods. Here we introduce WaterKit, a new grid-based sampling method with explicit water molecules to calculate thermodynamic properties using the GIST method. Our results show that the discrete placement of water molecules is successful in reproducing the position of crystallographic waters with very high accuracy, as well as providing thermodynamic estimates with accuracy comparable to more expensive MD simulations. Unlike these methods, WaterKit can be used to analyze specific regions on the protein surface, (such as the binding site of a receptor), without having to hydrate and simulate the whole receptor structure. The results show the feasibility of a general and fast method to compute thermodynamic properties of water molecules, making it well-suited to be integrated in high-throughput pipelines such as molecular docking.

Unveiling the Inhibitory Potentials of Peptidomimetic Azanitriles and Pyridyl Esters towards SARS-CoV-2 Main Protease: A Molecular Modelling Investigation

The severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) is responsible for COVID-19, which was declared a global pandemic in March 2020 by the World Health Organization (WHO). Since SARS-CoV-2 main protease plays an essential role in the virus’s life cycle, the design of small drug molecules with lower molecular weight has been a promising development targeting its inhibition. Herein, we evaluated the novel peptidomimetic azatripeptide and azatetrapeptide nitriles against SARS-CoV-2 main protease. We employed molecular dynamics (MD) simulations to elucidate the selected compounds’ binding free energy profiles against SARS-CoV-2 and further unveil the residues responsible for the drug-binding properties. Compound 8 exhibited the highest binding free energy of −49.37 ± 0.15 kcal/mol, followed by compound 7 (−39.83 ± 0.19 kcal/mol), while compound 17 showed the lowest binding free energy (−23.54 ± 0.19 kcal/mol). In addition, the absorption, distribution, metabolism, and excretion (ADME) assessment was performed and revealed that only compound 17 met the drug-likeness parameters and exhibited high pharmacokinetics to inhibit CYP1A2, CYP2C19, and CYP2C9 with better absorption potential and blood-brain barrier permeability (BBB) index. The additional intermolecular evaluations suggested compound 8 as a promising drug candidate for inhibiting SARS-CoV-2 Mpro. The substitution of isopropane in compound 7 with an aromatic benzene ring in compound 8 significantly enhanced the drug’s ability to bind better at the active site of the SARS-CoV-2 Mpro.

Enhanced Grand Canonical Sampling of Occluded Water Sites Using Nonequilibrium Candidate Monte Carlo

Water molecules play a key role in many biomolecular systems, particularly when bound at protein-ligand interfaces. However, molecular simulation studies on such systems are hampered by the relatively long time scales over which water exchange between a protein and solvent takes place. Grand canonical Monte Carlo (GCMC) is a simulation technique that avoids this issue by attempting the insertion and deletion of water molecules within a given structure. The approach is constrained by low acceptance probabilities for insertions in congested systems, however. To address this issue, here, we combine GCMC with nonequilibium candidate Monte Carlo (NCMC) to yield a method that we refer to as grand canonical nonequilibrium candidate Monte Carlo (GCNCMC), in which the water insertions and deletions are carried out in a gradual, nonequilibrium fashion. We validate this new approach by comparing GCNCMC and GCMC simulations of bulk water and three protein binding sites. We find that not only is the efficiency of the water sampling improved by GCNCMC but that it also results in increased sampling of ligand conformations in a protein binding site, revealing new water-mediated ligand-binding geometries that are not observed using alternative enhanced sampling techniques.

The slow but steady rise of binding free energy calculations in drug discovery

Binding free energy calculations are increasingly used in drug discovery research to predict protein-ligand binding affinities and to prioritize candidate drug molecules accordingly. It has taken decades of collective effort to transform this academic concept into a technology adopted by the pharmaceutical and biotech industry. Having personally witnessed and taken part in this transformation, here I recount the (incomplete) list of problems that had to be solved to make this computational tool practical and suggest areas of future development.

Unraveling Complex Local Protein Environments with 4-Cyano-l-phenylalanine

We present a multifaceted approach to effectively probe complex local protein environments utilizing the vibrational reporter unnatural amino acid (UAA) 4-cyano-l-phenylalanine (pCNPhe) in the model system superfolder green fluorescent protein (sfGFP). This approach combines temperature-dependent infrared (IR) spectroscopy, X-ray crystallography, and molecular dynamics (MD) simulations to provide a molecular interpretation of the local environment of the nitrile group in the protein. Specifically, a two-step enantioselective synthesis was developed that provided an 87% overall yield of pCNPhe in high purity without the need for chromatography. It was then genetically incorporated individually at three unique sites (74, 133, and 149) in sfGFP to probe these local protein environments. The incorporation of the UAA site-specifically in sfGFP utilized an engineered, orthogonal tRNA synthetase in E. coli using the Amber codon suppression protocol, and the resulting UAA-containing sfGFP constructs were then explored with this approach. This methodology was effectively utilized to further probe the local environments of two surface sites (sites 133 and 149) that we previously explored with room temperature IR spectroscopy and X-ray crystallography and a new interior site (site 74) featuring a complex local environment around the nitrile group of pCNPhe. Site 133 was found to be solvent-exposed, while site 149 was partially buried. Site 74 was found to consist of three distinct local environments around the nitrile group including nonspecific van der Waals interactions, hydrogen-bonding to a structural water, and hydrogen-bonding to a histidine side chain.

Forces Driving a Magic Bullet to Its Target: Revisiting the Role of Thermodynamics in Drug Design, Development, and Optimization

Drug discovery strategies have advanced significantly towards prioritizing target selectivity to achieve the longstanding goal of identifying "magic bullets" amongst thousands of chemical molecules screened for therapeutic efficacy. A myriad of emerging and existing health threats, including the SARS-CoV-2 pandemic, alarming increase in bacterial resistance, and potentially fatal chronic ailments, such as cancer, cardiovascular disease, and neurodegeneration, have incentivized the discovery of novel therapeutics in treatment regimens. The design, development, and optimization of lead compounds represent an arduous and time-consuming process that necessitates the assessment of specific criteria and metrics derived via multidisciplinary approaches incorporating functional, structural, and energetic properties. The present review focuses on specific methodologies and technologies aimed at advancing drug development with particular emphasis on the role of thermodynamics in elucidating the underlying forces governing ligand-target interaction selectivity and specificity. In the pursuit of novel therapeutics, isothermal titration calorimetry (ITC) has been utilized extensively over the past two decades to bolster drug discovery efforts, yielding information-rich thermodynamic binding signatures. A wealth of studies recognizes the need for mining thermodynamic databases to critically examine and evaluate prospective drug candidates on the basis of available metrics. The ultimate power and utility of thermodynamics within drug discovery strategies reside in the characterization and comparison of intrinsic binding signatures that facilitate the elucidation of structural-energetic correlations which assist in lead compound identification and optimization to improve overall therapeutic efficacy.

Absolute binding free energy calculations improve enrichment of actives in virtual compound screening

We determined the effectiveness of absolute binding free energy (ABFE) calculations to refine the selection of active compounds in virtual compound screening, a setting where the more commonly used relative binding free energy approach is not readily applicable. To do this, we conducted baseline docking calculations of structurally diverse compounds in the DUD-E database for three targets, BACE1, CDK2 and thrombin, followed by ABFE calculations for compounds with high docking scores. The docking calculations alone achieved solid enrichment of active compounds over decoys. Encouragingly, the ABFE calculations then improved on this baseline. Analysis of the results emphasizes the importance of establishing high quality ligand poses as starting points for ABFE calculations, a nontrivial goal when processing a library of diverse compounds without informative co-crystal structures. Overall, our results suggest that ABFE calculations can play a valuable role in the drug discovery process.

In Silico Drug Repurposing of FDA-Approved Drugs Highlighting Promacta as a Potential Inhibitor of H7N9 Influenza Virus

Influenza virus infections continue to be a significant and recurrent public health problem. Although vaccine efficacy varies, regular immunisation is the most effective method for suppressing the influenza virus. Antiviral drugs are available for influenza, although two of the four FDA-approved antiviral treatments have resulted in significant drug resistance. Therefore, new treatments are being sought to reduce the burden of flu-related illness. The time-consuming development of treatments for new and re-emerging diseases such as influenza and the high failure rate are increasing concerns. In this context, we used an in silico-based drug repurposing method to repurpose FDA-approved drugs as potential therapies against the H7N9 virus. To find potential inhibitors, a total of 2568 drugs were screened. Promacta, tucatinib, and lurasidone were identified as promising hits in the DrugBank database. According to the calculations of MM-GBSA, tucatinib (−54.11 kcal/mol) and Promacta (−56.20 kcal/mol) occupied the active site of neuraminidase with a higher binding affinity than the standard drug peramivir (−49.09 kcal/mol). Molecular dynamics (MD) simulation studies showed that the C-α atom backbones of the complexes of tucatinib and Promacta neuraminidase were stable throughout the simulation period. According to ADME analysis, the hit compounds have a high gastrointestinal absorption (GI) and do not exhibit properties that allow them to cross the blood–brain barrier (BBB). According to the in silico toxicity prediction, Promacta is not cardiotoxic, while lurasidone and tucatinib show only weak inhibition. Therefore, we propose to test these compounds experimentally against the influenza H7N9 virus. The investigation and validation of these potential H7N9 inhibitors would be beneficial in order to bring these compounds into clinical settings.

Structural Investigations and Binding Mechanisms of Oseltamivir Drug Resistance Conferred by the E119V Mutation in Influenza H7N9 Virus

The use of vaccinations and antiviral medications have gained popularity in the therapeutic management of avian influenza H7N9 virus lately. Antiviral medicines are more popular due to being readily available. The presence of the neuraminidase protein in the avian influenza H7N9 virus and its critical role in the cleavage of sialic acid have made it a target drug in the development of influenza virus drugs. Generally, the neuraminidase proteins have common conserved amino acid residues and any mutation that occurs around or within these conserved residues affects the susceptibility and replicability of the influenza H7N9 virus. Herein, we investigated the interatomic and intermolecular dynamic impacts of the experimentally reported E119V mutation on the oseltamivir resistance of the influenza H7N9 virus. We extensively employed molecular dynamic (MD) simulations and subsequent post-MD analyses to investigate the binding mechanisms of oseltamivir-neuraminidase wildtype and E119V mutant complexes. The results revealed that the oseltamivir-wildtype complex was more thermodynamically stable than the oseltamivir-E119V mutant complex. Oseltamivir exhibited a greater binding affinity for wildtype (−15.46 ± 0.23 kcal/mol) relative to the E119V mutant (−11.72 ± 0.21 kcal/mol). The decrease in binding affinity (−3.74 kcal/mol) was consistent with RMSD, RMSF, SASA, PCA, and hydrogen bonding profiles, confirming that the E119V mutation conferred lower conformational stability and weaker protein–ligand interactions. The findings of this oseltamivir-E119V mutation may further assist in the design of compounds to overcome E119V mutation in the treatment of influenza H7N9 virus patients.

Electronic Polarization Is Essential for the Stabilization and Dynamics of Buried Ion Pairs in Staphylococcal Nuclease Mutants

Buried charged residues play important roles in the modulation of protein stabilities and conformational dynamics and make crucial contributions to protein functions. Considering the generally nonpolar nature of protein interior, a key question concerns the contribution of electronic polarization to the stabilization and properties of buried charges. We answer this question by conducting free energy simulations using the latest polarizable CHARMM force field based on Drude oscillators for a series of Staphylococcal nuclease mutants that involve a buried Glu-Lys pair in different titration states and orientations. While a nonpolarizable model suggests that the ionized form of the buried Glu-Lys pair is more than 40 kcal/mol less stable than the charge-neutral form, the two titration states are comparable in stability when electronic polarization is included explicitly, a result better reconcilable with available experimental data. Analysis of free energy components suggests that additional stabilization of the ionized Glu-Lys pair has contributions from both the enhanced salt-bridge strength and stronger interaction between the ion-pair and surrounding protein residues and penetrated water. Despite the stronger direct interaction between Glu and Lys, the ion-pair exhibits considerably larger and faster structural fluctuations when polarization is included, due to compensation of interactions in the cavity. Collectively, observations from this work provide compelling evidence that electronic polarization is essential to the stability, hydration, dynamics, and therefore function of buried charges in proteins. Therefore, our study advocates for the explicit consideration of electronic polarization for mechanistic and engineering studies that implicate buried charged residues, such as enzymes and ion transporters.

Enhancing Sampling of Water Rehydration on Ligand Binding: A Comparison of Techniques

Water often plays a key role in protein structure, molecular recognition, and mediating protein-ligand interactions. Thus, free energy calculations must adequately sample water motions, which often proves challenging in typical MD simulation time scales. Thus, the accuracy of methods relying on MD simulations ends up limited by slow water sampling. Particularly, as a ligand is removed or modified, bulk water may not have time to fill or rearrange in the binding site. In this work, we focus on several molecular dynamics (MD) simulation-based methods attempting to help rehydrate buried water sites: BLUES, using nonequilibrium candidate Monte Carlo (NCMC); grand, using grand canonical Monte Carlo (GCMC); and normal MD. We assess the accuracy and efficiency of these methods in rehydrating target water sites. We selected a range of systems with varying numbers of waters in the binding site, as well as those where water occupancy is coupled to the identity or binding mode of the ligand. We analyzed the rehydration of buried water sites in binding pockets using both clustering of trajectories and direct analysis of electron density maps. Our results suggest both BLUES and grand enhance water sampling relative to normal MD and grand is more robust than BLUES, but also that water sampling remains a major challenge for all of the methods tested. The lessons we learned for these methods and systems are discussed.

Intermolecular Mechanism and Dynamic Investigation of Avian Influenza H7N9 Virus' Susceptibility to E119V-Substituted Peramivir-Neuraminidase Complex

The H7N9 virus attaches itself to the human cell receptor protein containing the polysaccharide that terminates with sialic acid. The mutation of neuraminidase at residue E119 has been explored experimentally. However, there is no adequate information on the substitution with E119V in peramivir at the intermolecular level. Therefore, a good knowledge of the interatomic interactions is a prerequisite in understanding its transmission mode and subsequent effective inhibitions of the sialic acid receptor cleavage by neuraminidase. Herein, we investigated the mechanism and dynamism on the susceptibility of the E119V mutation on the peramivir-neuraminidase complex relative to the wildtype complex at the intermolecular level. This study aims to investigate the impact of the 119V substitution on the neuraminidase-peramivir complex and unveil the residues responsible for the complex conformations. We employed molecular dynamic (MD) simulations and extensive post-MD analyses in the study. These extensive computational investigations were carried out on the wildtype and the E119V mutant complex of the protein for holistic insights in unveiling the effects of this mutation on the binding affinity and the conformational terrain of peramivir-neuraminidase E119V mutation. The calculated total binding energy (ΔGbind) for the peramivir wildtype is -49.09 ± 0.13 kcal/mol, while the E119V mutant is -58.55 ± 0.15 kcal/mol. The increase in binding energy (9.46 kcal/mol) is consistent with other post-MD analyses results, confirming that E119V substitution confers a higher degree of stability on the protein complex. This study promises to proffer contributory insight and additional knowledge that would enhance future drug designs and help in the fight targeted at controlling the avian influenza H7N9 virus. Therefore, we suggest that experimentalists collaborate with computational chemists for all investigations of this topic, as we have done in our previous studies.