Single-Molecule Pulling Simulations Can Discern Active from Inactive Enzyme Inhibitors

Examining sialic acid derivatives as potential inhibitors of SARS-CoV-2 spike protein receptor binding domain

Severe acute respiratory syndrome coronavirus 2 (SARS CoV-2) has been the primary reason behind the COVID-19 global pandemic which has affected millions of lives worldwide. The fundamental cause of the infection is the molecular binding of the viral spike protein receptor binding domain (SP-RBD) with the human cell angiotensin-converting enzyme 2 (ACE2) receptor. The infection can be prevented if the binding of RBD-ACE2 is resisted by utilizing certain inhibitors or drugs that demonstrate strong binding affinity towards the SP RBD. Sialic acid based glycans found widely in human cells and tissues have notable propensity of binding to viral proteins of the coronaviridae family. Recent experimental literature have used N-acetyl neuraminic acid (Sialic acid) to create diagnostic sensors for SARS-CoV-2, but a detailed interrogation of the underlying molecular mechanisms is warranted. Here, we perform all atom molecular dynamics (MD) simulations for the complexes of certain Sialic acid-based molecules with that of SP RBD of SARS CoV-2. Our results indicate that Sialic acid not only reproduces a binding affinity comparable to the RBD-ACE2 interactions, it also assumes the longest time to dissociate completely from the protein binding pocket of SP RBD. Our predictions corroborate that a combination of electrostatic and van der Waals energies as well the polar hydrogen bond interactions between the RBD residues and the inhibitors influence free energy of binding.Communicated by Ramaswamy H. Sarma.

Simulation of the ligand-leaving process of the human heat shock protein

The human heat shock protein plays a critical role in various diseases and is an important target for pharmacological modulation. Simulation of conformational changes and free energy profiles of the human heat shock protein derived by the ligand-leaving process is a challenging issue. In this work, steered molecular dynamics simulation was adopted to simulate the ligand-leaving process. Two composite systems of heat shock protein NHSP90 and small molecules 6FJ and 6G7 are selected as research objects. The free energy during the leaving of ligand small molecules is calculated using conventional molecular dynamics simulation, steered molecular dynamics simulation (SMD), and the umbrella sampling method. We found that the a slower pulling velocity (0.001 nm ns-1) will result in 2.19 kcal mol-1, and the umbrella sampling method gives a value of 3.26 kcal mol-1 for the free energy difference for the two systems, which reasonably agrees with experimental results. A faster-pulling velocity (0.01 nm ns-1) leads to a large overestimation of free energy. At the same time, the conformational analysis indicated that the faster pulling velocity may lead to the conformational change of NHSP90, which was proved to be false by the slower pulling velocity and the umbrella sampling method.

Loop Evolutionary Patterns Shape Catalytic Efficiency of TRI101/201 for Trichothecenes: Insights into Protein-Substrate Interactions

Trichothecenes are highly toxic mycotoxins produced by Fusarium fungi, while TRI101/201 family enzymes play a crucial role in detoxification through acetylation. Studies on the substrate specificity and catalytic kinetics of TRI101/201 have revealed distinct kinetic characteristics, with significant differences observed in catalytic efficiency toward deoxynivalenol, while the catalytic efficiency for T-2 toxin remains relatively consistent. In this study, we used structural bioinformatics analysis and a molecular dynamics simulation workflow to investigate the mechanism underlying the differential catalytic activity of TRI101/201. The findings revealed that the binding stability between trichothecenes and TRI101/201 hinges primarily on a hydrophobic cage structure within the binding site. An intrinsic disordered loop, termed loop cover, defined the evolutionary patterns of the TRI101/201 protein family that are categorized into four subfamilies (V1/V2/V3/M). Furthermore, the unique loop displayed different conformations among these subfamilies' structures, which served to disrupt (V1/V2/V3) or reinforce (M) the hydrophobic cages. The disrupted cages enhanced the water exposure of the hydrophilic moieties of substrates like deoxynivalenol and thereby hindered their binding to the catalytic sites of V-type enzymes. In contrast, this water exposure does not affect substrates like T-2 toxin, which have more hydrophobic substituents, resulting in a comparable catalytic efficiency of both V- and M-type enzymes. Overall, our studies provide theoretical support for understanding the catalytic mechanism of TRI101/201, which shows how an intrinsic disordered loop could impact the protein-ligand binding and suggests a direction for rational protein design in the future.

SARS-CoV-2 external structures interacting with nanospheres using docking and molecular dynamics

Coronavirus is caused by the SARS-CoV-2 virus has shown rapid proliferation and scarcity of treatments with proven effectiveness. In this way, we simulated the hospitalization of carbon nanospheres, with external active sites of the SARS-CoV-2 virus (M-Pro, S-Gly and E-Pro), which can be adsorbed or inactivated when interacting with the nanospheres. The computational procedures performed in this work were developed with the SwissDock server for molecular docking and the GROMACS software for molecular dynamics, making it possible to extract relevant data on affinity energy, distance between molecules, free Gibbs energy and mean square deviation of atomic positions, surface area accessible to solvents. Molecular docking indicates that all ligands have an affinity for the receptor's active sites. The nanospheres interact favorably with all proteins, showing promising results, especially C60, which presented the best affinity energy and RMSD values ​​for all protein macromolecules investigated. The C60 with E-Pro exhibited the highest affinity energy of -9.361 kcal/mol, demonstrating stability in both molecular docking and molecular dynamics simulations. Our RMSD calculations indicated that the nanospheres remained predominantly stable, fluctuating within a range of 2 to 3 Å. Additionally, the analysis of other structures yielded promising results that hold potential for application in other proteases.Communicated by Ramaswamy H. Sarma.

15 Years of molecular simulation of drug-binding kinetics

INTRODUCTION

Drug-binding kinetics has been increasingly recognized as an important factor to be considered in drug discovery. Long residence time could prolong the action of some drugs while produce toxicity on others. Early evaluation of the binding kinetics of drug candidates could reduce attrition rate late in the drug discovery process. Computational prediction of drug-binding kinetics is useful as compounds can be evaluated even before they are made. However, simulation of drug-binding kinetics is a challenging problem because of the long-time scale involved. Nevertheless, significant progress has been made.

AREAS COVERED

This review illustrates the rapid evolution of qualitative to quantitative molecular dynamics-based methods that have been developed over the last 15 years.

EXPERT OPINION

The development of new methods based on molecular dynamics simulations now enables computation of absolute association/dissociation rate constants. Cheaper methods capable of identifying candidates with fast or slow binding kinetics, or rank-ordering rate constants are also available. Together, these methods have generated useful insights into the molecular mechanisms of drug-binding kinetics, and the design of drug candidates with therapeutically favorable kinetics. Although predicting absolute rate constants is still expensive and challenging, rapid improvement is expected in the coming years with the continuing refinement of current technologies, development of new methodologies, and the utilization of machine learning.

Lenalidomide Stabilizes Protein-Protein Complexes by Turning Labile Intermolecular H-Bonds into Robust Interactions

Targeted protein degradation is a promising therapeutic strategy, spearheaded by the anti-myeloma drugs lenalidomide and pomalidomide. These drugs stabilize very efficiently the complex between the E3 ligase Cereblon (CRBN) and several non-native client proteins (neo-substrates), including the transcription factors Ikaros and Aiolos and the enzyme Caseine Kinase 1α (CK1α,), resulting in their degradation. Although the structures for these complexes have been determined, there are no evident interactions that can account for the high efficiency of formation of the ternary complex. We show that lenalidomide's stabilization of the CRBN-CK1α complex is largely due to hydrophobic shielding of intermolecular hydrogen bonds. We also find a quantitative relationship between hydrogen bond robustness and binding affinities of the ternary complexes. These results pave the way to further understand cooperativity effects in drug-induced protein-protein complexes and could help in the design of improved molecular glues and more efficient protein degraders.

Korupensamine A, but not its atropisomer, korupensamine B, inhibits SARS-CoV-2 in vitro by targeting its main protease (Mpro)

By combining docking and molecular dynamics simulations, we explored a library of 65 mostly axially chiral naphthylisoquinoline alkaloids and their analogues, with most different molecular architectures and structural analogues, for their activity against SARS-CoV-2. Although natural biaryls are often regarded without consideration of their axial chirality, they can bind to protein targets in an atroposelective manner. By combining docking results with steered molecular dynamics simulations, we identified one alkaloid, korupensamine A, that atropisomer-specifically inhibited the main protease (M^pro) activity of SARS-CoV-2 significantly in comparison to the reference covalent inhibitor GC376 (IC₅₀ = 2.52 ± 0.14 and 0.88 ± 0.15 μM, respectively) and reduced viral growth by five orders of magnitude in vitro (EC₅₀ = 4.23 ± 1.31 μM). To investigate the binding pathway and mode of interaction of korupensamine A within the active site of the protease, we utilized Gaussian accelerated molecular dynamics simulations, which reproduced the docking pose of korupensamine A inside the active site of the enzyme. The study presents naphthylisoquinoline alkaloids as a new class of potential anti-COVID-19 agents.

Implementation of Telescoping Boxes in Adaptive Steered Molecular Dynamics

Long-time dynamical processes, such as those involving protein unfolding and ligand interactions, can be accelerated and realized through steered molecular dynamics (SMD). The challenge has been the extraction of information from such simulations that generalize for complex nonequilibrium processes. The use of Jarzynski's equality opened the possibility of determining the free energy along the steered coordinate, but sampling over the nonequilibrium trajectories is slow to converge. Adaptive steered molecular dynamics (ASMD) and other related techniques have been introduced to overcome this challenge through the use of stages. Here, we take advantage of these stages to address the numerical cost that arises from the required use of very large solvent boxes. We introduce telescoping box schemes within adaptive steered molecular dynamics (ASMD) in which we adjust the solvent box between stages and thereby vary (and optimize) the required number of solvent molecules. We have benchmarked the method on a relatively long α-helical peptide, Ala₃₀, with respect to the potential of mean force and hydrogen bonds. We show that the use of telescoping boxes introduces little numerical error while significantly reducing the computational cost.

Determination of Multidirectional Pathways for Ligand Release from the Receptor: A New Approach Based on Differential Evolution

Steered molecular dynamics (SMD) simulation is a powerful method in computer-aided drug design as it can be used to access the relative binding affinity with high precision but with low computational cost. The success of SMD depends on the choice of the direction along which the ligand is pulled from the receptor-binding site. In most simulations, the unidirectional pathway was used, but in some cases, this choice resulted in the ligand colliding with the complex surface of the exit tunnel. To overcome this difficulty, several variants of SMD with multidirectional pulling have been proposed, but they are not completely devoid of disadvantages. Here, we have proposed to determine the direction of pulling with a simple scoring function that minimizes the receptor–ligand interaction, and an optimization algorithm called differential evolution is used for energy minimization. The effectiveness of our protocol was demonstrated by finding expulsion pathways of Huperzine A and camphor from the binding site of Torpedo California acetylcholinesterase and P450cam proteins, respectively, and comparing them with the previous results obtained using memetic sampling and random acceleration molecular dynamics. In addition, by applying this protocol to a set of ligands bound with LSD1 (lysine specific demethylase 1), we obtained a much higher correlation between the work of pulling force and experimental data on the inhibition constant IC50 compared to that obtained using the unidirectional approach based on minimal steric hindrance.

Impact of PROTAC Linker Plasticity on the Solution Conformations and Dissociation of the Ternary Complex

The conformational behavior of a small molecule free in solution is important to understand the free energy of binding to its target. This could be of special interest for proteolysis-targeting chimeras (PROTACs) due to their often flexible and lengthy linkers and the need to induce a ternary complex. Here, we report on the molecular dynamics (MD) simulations of two PROTACs, MZ1 and dBET6, revealing different linker conformational behaviors. The simulation of MZ1 in dimethyl sulfoxide (DMSO) agrees well with the nuclear magnetic resonance study, providing strong support for the relevance of our simulations. To further understand the role of linker plasticity in the formation of a ternary complex, the dissociation of the complex von Hippel-Lindau-MZ1-BRD4 is investigated in detail by steered simulations and is shown to follow a two-step pathway. Interestingly, both MZ1 and dBET6 display in water, a tendency toward an intramolecular lipophilic interaction between the two warheads. The hydrophobic contact of the two warheads would prevent them from binding to their respective proteins and might have an effect on the efficacy of induced cellular protein degradation. However, conformations featuring this hydrophobic contact of the two warheads are calculated to be marginally more favorable.

The prediction of protein-ligand unbinding for modern drug discovery

INTRODUCTION

Drug-target thermodynamic and kinetic information have perennially important roles in drug design. The prediction of protein-ligand unbinding, which can provide important kinetic information, in experiments continues to face great challenges. Uncovering protein-ligand unbinding through molecular dynamics simulations has become efficient and inexpensive with the progress and enhancement of computing power and sampling methods.

AREAS COVERED

In this review, various sampling methods for protein-ligand unbinding and their basic principles are firstly briefly introduced. Then, their applications in predicting aspects of protein-ligand unbinding, including unbinding pathways, dissociation rate constants, residence time and binding affinity, are discussed.

EXPERT OPINION

Although various sampling methods have been successfully applied in numerous systems, they still have shortcomings and deficiencies. Most enhanced sampling methods require researchers to possess a wealth of prior knowledge of collective variables or reaction coordinates. In addition, most systems studied at present are relatively simple, and the study of complex systems in real drug research remains greatly challenging. Through the combination of machine learning and enhanced sampling methods, prediction accuracy can be further improved, and some problems encountered in complex systems also may be solved.

Dependence of Work on the Pulling Speed in Mechanical Ligand Unbinding

In single-molecule force spectroscopy, the rupture force F_max required for mechanical unfolding of a biomolecule or for pulling a ligand out of a binding site depends on the pulling speed V and, in the linear Bell-Evans regime, F_max ∼ ln(V). Recently, it has been found that non-equilibrium work W is better than F_max in describing relative ligand binding affinity, but the dependence of W on V remains unknown. In this paper, we developed an analytical theory showing that in the linear regime, W ∼ c₁ ln(V) + c₂ ln²(V), where c₁ and c₂ are constants. This quadratic dependence was also confirmed by all-atom steered molecular dynamics simulations of protein-ligand complexes. Although our theory was developed for ligand unbinding, it is also applicable to other processes, such as mechanical unfolding of proteins and other biomolecules, due to its universality.

Metadynamics-Based Approaches for Modeling the Hypoxia-Inducible Factor 2α Ligand Binding Process

Several methods based on enhanced-sampling molecular dynamics have been proposed for studying ligand binding processes. Here, we developed a protocol that combines the advantages of steered molecular dynamics (SMD) and metadynamics. While SMD is proposed for investigating possible unbinding pathways of the ligand and identifying the preferred one, metadynamics, with the path collective variable (PCV) formalism, is suggested to explore the binding processes along the pathway defined on the basis of SMD, by using only two CVs. We applied our approach to the study of binding of two known ligands to the hypoxia-inducible factor 2α, where the buried binding cavity makes simulation of the process a challenging task. Our approach allowed identification of the preferred entrance pathway for each ligand, highlighted the features of the bound and intermediate states in the free-energy surface, and provided a binding affinity scale in agreement with experimental data. Therefore, it seems to be a suitable tool for elucidating ligand binding processes of similar complex systems.

Discovery of RTA ricin subunit inhibitors: a computational study using PM7 quantum chemical method and steered molecular dynamics

Ricin is a potent toxin derived from the castor bean plant and comprises two subunits, RTA and RTB. Because of its cytotoxicity, ricin has alarmed world authorities for its potential use as a chemical weapon. Ricin also affects castor bean agribusiness, given the risk of animal and human poisoning. Over the years, many groups attempted to propose small-molecules that bind to the RTA active site, the catalytic chain. Despite such efforts, there is still no effective countermeasure against ricin poisoning. The computational study carried out in the present work renews the discussion about small-molecules that may inhibit this toxin. Here, a structure-based virtual screening protocol capable of discerning active RTA inhibitors from inactive ones was performed to screen over 2 million compounds from the ZINC database to find novel scaffolds that strongly bind into the active site of the RTA. Besides, a novel score method based on ligand undocking force profiles and semi-empirical quantum chemical calculations provided insights into the rescore of docking poses. Summing up, the filtering steps pointed out seven main compounds, with the SCF00-451 as a promising candidate to inhibit the killing activity of such potent phytotoxin.

Thermodynamics and Kinetics of Drug-Target Binding by Molecular Simulation

Computational studies play an increasingly important role in chemistry and biophysics, mainly thanks to improvements in hardware and algorithms. In drug discovery and development, computational studies can reduce the costs and risks of bringing a new medicine to market. Computational simulations are mainly used to optimize promising new compounds by estimating their binding affinity to proteins. This is challenging due to the complexity of the simulated system. To assess the present and future value of simulation for drug discovery, we review key applications of advanced methods for sampling complex free-energy landscapes at near nonergodicity conditions and for estimating the rate coefficients of very slow processes of pharmacological interest. We outline the statistical mechanics and computational background behind this research, including methods such as steered molecular dynamics and metadynamics. We review recent applications to pharmacology and drug discovery and discuss possible guidelines for the practitioner. Recent trends in machine learning are also briefly discussed. Thanks to the rapid development of methods for characterizing and quantifying rare events, simulation's role in drug discovery is likely to expand, making it a valuable complement to experimental and clinical approaches.

Quantifying entropic barriers in single-molecule surface diffusion

The quantitative role of entropy in the surface diffusion of molecules with many degrees of freedom is still not well understood. Here, we quantify entropic diffusion barriers as well as attempt frequencies by performing a systematic decomposition of the Arrhenius equation for single oligophenyl molecules of various lengths (two to six phenyl rings and benzene as the reference) on an amorphous silica surface using extensive molecular dynamics simulations. Attempt frequencies evaluated from velocity auto-correlation functions are found close to k_BT/h, the frequency factor of transition state theory. Importantly, we find large positive entropy contributions to the free energy barrier of diffusion up to 55%, increasing with molecular length with 4.1 kJ/mol/phenyl ring. The entropic barrier is about 40%-60% of the entropy of the molecule surface adsorption free energy, revealing that at the transition states, the molecules can liberate a major part of their conformational states, increasing with length. The substantial role of the internal degrees of freedom for the diffusive dynamics is explicitly demonstrated by studying internally constrained, "rigid" version of the molecules. Finally, we discuss also rotational diffusion and the role of surface vibrations. Our results affirm that it is essential for quantitative studies and interpretation of surface diffusion of complex molecules to consider internal entropic effects.

Merging Ligand-Based and Structure-Based Methods in Drug Discovery: An Overview of Combined Virtual Screening Approaches

Virtual screening (VS) is an outstanding cornerstone in the drug discovery pipeline. A variety of computational approaches, which are generally classified as ligand-based (LB) and structure-based (SB) techniques, exploit key structural and physicochemical properties of ligands and targets to enable the screening of virtual libraries in the search of active compounds. Though LB and SB methods have found widespread application in the discovery of novel drug-like candidates, their complementary natures have stimulated continued efforts toward the development of hybrid strategies that combine LB and SB techniques, integrating them in a holistic computational framework that exploits the available information of both ligand and target to enhance the success of drug discovery projects. In this review, we analyze the main strategies and concepts that have emerged in the last years for defining hybrid LB + SB computational schemes in VS studies. Particularly, attention is focused on the combination of molecular similarity and docking, illustrating them with selected applications taken from the literature.

Metadynamics as a Postprocessing Method for Virtual Screening with Application to the Pseudokinase Domain of JAK2

With standard scoring methods, top-ranked compounds from virtual screening by docking often turn out to be inactive. For this reason, metadynamics, a method used to sample rare events, was studied to further evaluate docking poses with the aim of reducing false positives. Specifically, virtual screening was performed with Glide SP to seek potential molecules to bind to the ATP site in the pseudokinase domain of JAK2 kinase, and promising compounds were selected from the top-ranked 1000 based on visualization. Rescoring with Glide XP, GOLD, and MM/GBSA was unable to differentiate well between active and inactive compounds. Metadynamics was then used to gauge the relative binding affinity from the required time or the potential of mean force needed to dissociate the ligand from the bound complex. With consideration of previously known binders of varying affinities, metadynamics was able to differentiate between the most active compounds and inactive or weakly active ones, and it could identify correctly most of the selected virtual screening compounds as false positives. Thus, metadynamics has the potential to be a viable postprocessing method for virtual screening, minimizing the expense of buying or synthesizing inactive compounds.

MovableType Software for Fast Free Energy-Based Virtual Screening: Protocol Development, Deployment, Validation, and Assessment

For decades, the complicated energy surfaces found in macromolecular protein:ligand structures, which require large amounts of computational time and resources for energy state sampling, have been an inherent obstacle to fast, routine free energy estimation in industrial drug discovery efforts. Beginning in 2013, the Merz research group addressed this cost with the introduction of a novel sampling methodology termed “Movable Type” (MT). Using numerical integration methods, the MT method reduces the computational expense for energy state sampling by independently calculating each atomic partition function from an initial molecular conformation in order to estimate the molecular free energy using ensembles of the atomic partition functions. In this work, we report a software package, the DivCon Discovery Suite with the MovableType module from QuantumBio Inc., that performs this MT free energy estimation protocol in a fast, fully encapsulated manner. We discuss the computational procedures and improvements to the original work, and we detail the corresponding settings for this software package. Finally, we introduce two validation benchmarks to evaluate the overall robustness of the method against a broad range of protein:ligand structural cases. With these publicly available benchmarks, we show that the method can use a variety of input types and parameters and exhibits comparable predictability whether the method is presented with “expensive” X-ray structures or “inexpensively docked” theoretical models. We also explore some next steps for the method. The MovableType software is available at http://www.quantumbioinc.com/

Identification of 14 Known Drugs as Inhibitors of the Main Protease of SARS-CoV-2

A consensus virtual screening protocol has been applied to ca. 2000 approved drugs to seek inhibitors of the main protease (M ^pro ) of SARS-CoV-2, the virus responsible for COVID-19. 42 drugs emerged as top candidates, and after visual analyses of the predicted structures of their complexes with M ^pro , 17 were chosen for evaluation in a kinetic assay for M ^pro inhibition. Remarkably 14 of the compounds at 100-μM concentration were found to reduce the enzymatic activity and 5 provided IC ₅₀ values below 40 μM: manidipine (4.8 μM), boceprevir (5.4 μM), lercanidipine (16.2 μM), bedaquiline (18.7 μM), and efonidipine (38.5 μM). Structural analyses reveal a common cloverleaf pattern for the binding of the active compounds to the P1, P1', and P2 pockets of M ^pro . Further study of the most active compounds in the context of COVID-19 therapy is warranted, while all of the active compounds may provide a foundation for lead optimization to deliver valuable chemotherapeutics to combat the pandemic.

Does SARS-CoV-2 Bind to Human ACE2 More Strongly Than Does SARS-CoV?

The 2019 novel coronavirus (SARS-CoV-2) epidemic, which was first reported in December 2019 in Wuhan, China, was declared a pandemic by the World Health Organization in March 2020. Genetically, SARS-CoV-2 is closely related to SARS-CoV, which caused a global epidemic with 8096 confirmed cases in more than 25 countries from 2002 to 2003. Given the significant morbidity and mortality rate, the current pandemic poses a danger to all of humanity, prompting us to understand the activity of SARS-CoV-2 at the atomic level. Experimental studies have revealed that spike proteins of both SARS-CoV-2 and SARS-CoV bind to angiotensin-converting enzyme 2 (ACE2) before entering the cell for replication. However, the binding affinities reported by different groups seem to contradict each other. Wrapp et al. (Science2020, 367, 1260–1263) showed that the spike protein of SARS-CoV-2 binds to the ACE2 peptidase domain (ACE2-PD) more strongly than does SARS-CoV, and this fact may be associated with a greater severity of the new virus. However, Walls et al. (Cell2020, 181, 281–292) reported that SARS-CoV-2 exhibits a higher binding affinity, but the difference between the two variants is relatively small. To understand the binding mechnism and experimental results, we investigated how the receptor binding domain (RBD) of SARS-CoV (SARS-CoV-RBD) and SARS-CoV-2 (SARS-CoV-2-RBD) interacts with a human ACE2-PD using molecular modeling. We applied a coarse-grained model to calculate the dissociation constant and found that SARS-CoV-2 displays a 2-fold higher binding affinity. Using steered all-atom molecular dynamics simulations, we demonstrate that, like a coarse-grained simulation, SARS-CoV-2-RBD was associated with ACE2-PD more strongly than was SARS-CoV-RBD, as evidenced by a higher rupture force and larger pulling work. We show that the binding affinity of both viruses to ACE2 is driven by electrostatic interactions.

Targeting the Protein Tunnels of the Urease Accessory Complex: A Theoretical Investigation

Urease is a nickel-containing enzyme that is essential for the survival of several and often deadly pathogenic bacterial strains, including Helicobacter pylori. Notwithstanding several attempts, the development of direct urease inhibitors without side effects for the human host remains, to date, elusive. The recently solved X-ray structure of the HpUreDFG accessory complex involved in the activation of urease opens new perspectives for structure-based drug discovery. In particular, the quaternary assembly and the presence of internal tunnels for nickel translocation offer an intriguing possibility to target the HpUreDFG complex in the search of indirect urease inhibitors. In this work, we adopted a theoretical framework to investigate such a hypothesis. Specifically, we searched for putative binding sites located at the protein–protein interfaces on the HpUreDFG complex, and we challenged their druggability through structure-based virtual screening. We show that, by virtue of the presence of tunnels, some protein–protein interfaces on the HpUreDFG complex are intrinsically well suited for hosting small molecules, and, as such, they possess good potential for future drug design endeavors.

How Good is Jarzynski's Equality for Computer-Aided Drug Design?

Accurate determination of the binding affinity of the ligand to the receptor remains a difficult problem in computer-aided drug design. Here, we study and compare the efficiency of Jarzynski’s equality (JE) combined with steered molecular dynamics and the linear interaction energy (LIE) method by assessing the binding affinity of 23 small compounds to six receptors, including β-lactamase, thrombin, factor Xa, HIV-1 protease (HIV), myeloid cell leukemia-1, and cyclin-dependent kinase 2 proteins. It was shown that Jarzynski’s nonequilibrium binding free energy ΔG_neq^Jar correlates with the available experimental data with the correlation levels R = 0.89, 0.86, 0.83, 0.80, 0.83, and 0.81 for six data sets, while for the binding free energy ΔG_LIE obtained by the LIE method, we have R = 0.73, 0.80, 0.42, 0.23, 0.85, and 0.01. Therefore, JE is recommended to be used for ranking binding affinities as it provides accurate and robust results. In contrast, LIE is not as reliable as JE, and it should be used with caution, especially when it comes to new systems.

Use of Molecular Dynamics Simulations in Structure-Based Drug Discovery

BACKGROUND

Traditional drug discovery is a lengthy process which involves a huge amount of resources. Modern-day drug discovers various multidisciplinary approaches amongst which, computational ligand and structure-based drug designing methods contribute significantly. Structure-based drug designing techniques require the knowledge of structural information of drug target and drug-target complexes. Proper understanding of drug-target binding requires the flexibility of both ligand and receptor to be incorporated. Molecular docking refers to the static picture of the drug-target complex(es). Molecular dynamics, on the other hand, introduces flexibility to understand the drug binding process.

OBJECTIVE

The aim of the present study is to provide a systematic review on the usage of molecular dynamics simulations to aid the process of structure-based drug design.

METHOD

This review discussed findings from various research articles and review papers on the use of molecular dynamics in drug discovery. All efforts highlight the practical grounds for which molecular dynamics simulations are used in drug designing program. In summary, various aspects of the use of molecular dynamics simulations that underline the basis of studying drug-target complexes were thoroughly explained.

RESULTS

This review is the result of reviewing more than a hundred papers. It summarizes various problems that use molecular dynamics simulations.

CONCLUSION

The findings of this review highlight how molecular dynamics simulations have been successfully implemented to study the structure-function details of specific drug-target complexes. It also identifies the key areas such as stability of drug-target complexes, ligand binding kinetics and identification of allosteric sites which have been elucidated using molecular dynamics simulations.

Enhanced sampling and free energy calculations for protein simulations

Molecular dynamics simulation is a powerful computational technique to study biomolecular systems, which complements experiments by providing insights into the structural dynamics relevant to biological functions at atomic scale. It can also be used to calculate the free energy landscapes of the conformational transitions to better understand the functions of the biomolecules. However, the sampling of biomolecular configurations is limited by the free energy barriers that need to be overcome, leading to considerable gaps between the timescales reached by MD simulation and those governing biological processes. To address this issue, many enhanced sampling methodologies have been developed to increase the sampling efficiency of molecular dynamics simulations and free energy calculations. Usually, enhanced sampling algorithms can be classified into methods based on collective variables (CV-based) and approaches which do not require predefined CVs (CV-free). In this chapter, the theoretical basis of free energy estimation is briefly reviewed first, followed by the reviews of the most common CV-based and CV-free methods including the presentation of some examples and recent developments. Finally, the combination of different enhanced sampling methods is discussed.

Rigorous Computational and Experimental Investigations on MDM2/MDMX-Targeted Linear and Macrocyclic Peptides

There is interest in peptide drug design, especially for targeting intracellular protein-protein interactions. Therefore, the experimental validation of a computational platform for enabling peptide drug design is of interest. Here, we describe our peptide drug design platform (CMDInventus) and demonstrate its use in modeling and predicting the structural and binding aspects of diverse peptides that interact with oncology targets MDM2/MDMX in comparison to both retrospective (pre-prediction) and prospective (post-prediction) data. In the retrospective study, CMDInventus modules (CMDpeptide, CMDboltzmann, CMDescore and CMDyscore) were used to accurately reproduce structural and binding data across multiple MDM2/MDMX data sets. In the prospective study, CMDescore, CMDyscore and CMDboltzmann were used to accurately predict binding affinities for an Ala-scan of the stapled α-helical peptide ATSP-7041. Remarkably, CMDboltzmann was used to accurately predict the results of a novel D-amino acid scan of ATSP-7041. Our investigations rigorously validate CMDInventus and support its utility for enabling peptide drug design.

Steering the Lipid Transfer To Unravel the Mechanism of Cholesteryl Ester Transfer Protein Inhibition

Human plasma cholesteryl ester transfer protein (CETP) mediates the transfer of neutral lipids from antiatherogenic high-density lipoproteins (HDLs) to proatherogenic low-density lipoproteins (LDLs). Recent cryo-electron microscopy studies have suggested that CETP penetrates its N- and C-terminal domains in HDL and LDL to form a ternary complex, which facilitates the lipid transfer between different lipoproteins. Inhibition of CETP lipid transfer activity has been shown to increase the plasma HDL-C levels and, therefore, became an effective strategy for combating cardiovascular diseases. Thus, understanding the molecular mechanism of inhibition of lipid transfer through CETP is of paramount importance. Recently reported inhibitors, torcetrapib and anacetrapib, exhibited low potency in addition to severe side effects, which essentially demanded a thorough knowledge of the inhibition mechanism. Here, we employ steered molecular dynamics simulations to understand how inhibitors interfere with the neutral lipid transfer mechanism of CETP. Our study revealed that inhibitors physically occlude the tunnel posing a high energy barrier for lipid transfer. In addition, inhibitors bring about the conformational changes in CETP that hamper CE passage and expose protein residues that disrupt the optimal hydrophobicity of the CE transfer path. The atomic level details presented here could accelerate the designing of safe and efficacious CETP inhibitors.

Unbinding of Translocator Protein 18 kDa (TSPO) Ligands: From in Vitro Residence Time to in Vivo Efficacy via in Silico Simulations

Translocator protein 18 kDa (TSPO) is a validated pharmacological target for the development of new treatments for neurological disorders. N,N-Dialkyl-2-phenylindol-3-ylglyoxylamides (PIGAs) are effective TSPO modulators and potentially useful therapeutics for the treatment of anxiety, central nervous system pathologies featuring astrocyte loss, and inflammatory-based neuropathologies. For this class of compounds, no correlation exists between the TSPO binding affinity and the corresponding functional efficacy. Rather, their biological effectiveness correlates with the kinetics of the unbinding events and more specifically with the residence time (RT). So far, the structural reasons for the different recorded RT of congeneric PIGAs remain elusive. Here, to understand the different kinetics of PIGAs, their unbinding paths were studied by employing enhanced-sampling molecular dynamics simulations. Results of these studies revealed how subtle structural differences between PIGAs have a substantial effect on the unbinding energetics. In particular, during the egress from the TSPO binding site, slow-dissociating PIGAs find tight interactions with the protein LP1 region thereby determining a long RT. Further support to these findings was achieved by in vivo studies, which demonstrated how the anxiolytic effect observed for the inspected PIGAs correlated with their RT to TSPO.

Computational Study for the Unbinding Routes of β-N-Acetyl-d-Hexosaminidase Inhibitor: Insight from Steered Molecular Dynamics Simulations

β-N-Acetyl-d-hexosaminidase from Ostrinia furnacalis (OfHex1) is a new target for the design of insecticides. Although some of its inhibitors have been found, there is still no commercial drug available at present. The residence time of the ligand may be important for its pharmacodynamic effect. However, the unbinding routes of ligands from OfHex1 still remain largely unexplored. In the present study, we first simulated the six dissociation routes of N,N,N-trimethyl-d-glucosamine-chitotriomycin (TMG-chitotriomycin, a highly selective inhibitor of OfHex1) from the active pocket of OfHex1 by steered molecular dynamics simulations. By comparing the potential of mean forces (PMFs) of six routes, Route 1 was considered as the most possible route with the lowest energy barrier. Furthermore, the structures of six different states for Route 1 were snapshotted, and the key amino acid residues affecting the dissociated time were analyzed in the unbinding pathway. Moreover, we also analyzed the "open⁻close" mechanism of Glu368 and Trp448 and found that their conformational changes directly affected the dissociation of TMG-chitotriomycin. Our findings would be helpful to understanding and identifying novel inhibitors against OfHex1 from virtual screening or lead-optimization.

Solvents to Fragments to Drugs: MD Applications in Drug Design

Simulations of molecular dynamics (MD) are playing an increasingly important role in structure-based drug discovery (SBDD). Here we review the use of MD for proteins in aqueous solvation, organic/aqueous mixed solvents (MDmix) and with small ligands, to the classic SBDD problems: Binding mode and binding free energy predictions. The simulation of proteins in their condensed state reveals solvent structures and preferential interaction sites (hot spots) on the protein surface. The information provided by water and its cosolvents can be used very effectively to understand protein ligand recognition and to improve the predictive capability of well-established methods such as molecular docking. The application of MD simulations to the study of the association of proteins with drug-like compounds is currently only possible for specific cases, as it remains computationally very expensive and labor intensive. MDmix simulations on the other hand, can be used systematically to address some of the common tasks in SBDD. With the advent of new tools and faster computers we expect to see an increase in the application of mixed solvent MD simulations to a plethora of protein targets to identify new drug candidates.

Probing the binding mode and unbinding mechanism of LSD1 inhibitors by combined computational methods

Lysine specific demethylase 1 (LSD1) has emerged as a potential drug target in cancer therapy and a variety of inhibitors have been reported. We have recently reported the discovery of a series of triazole-dithiocarbamate based compounds, which were basically confirmed as cofactor flavin adenine dinucleotide (FAD)-competing inhibitors by experiments. However, the binding modes of the inhibitors to the binding site were undetermined. Here, we employed computational methods including molecular docking, classical molecular dynamics (MD) and steered molecular dynamics (SMD) simulations to investigate the potential binding modes of these inhibitors to LSD1. Based on the high correlation between the mean non-equilibrium pulling work W and experimental binding affinity, we identified the optimal binding modes of this class of compounds with LSD1. Using the optimal inhibitor binding conformation, we then performed SMD to study the ligand unbinding mechanism with a lower pulling velocity at 0.0005 nm ps-1. We found that residue Arg316 plays a crucial role in the binding/unbinding process. Furthermore, a gatekeeper residue Trp756 influences the ligand unbinding process by acting like a switch via steric hindrance but can enhance the hydrophobic interaction with the inhibitor. Hydrophobic interaction also dominated the interaction between LSD1 and the inhibitors. The pivotal residues and interactions between LSD1 and inhibitors determined from this study can be used to improve the inhibition activity of this series of inhibitors in development and to discover new scaffolds as FAD-competing inhibitors in compound screening.

A Multiscale Simulation Approach to Modeling Drug-Protein Binding Kinetics

Drug-target binding kinetics has recently emerged as a sometimes critical determinant of in vivo efficacy and toxicity. Its rational optimization to improve potency or reduce side effects of drugs is, however, extremely difficult. Molecular simulations can play a crucial role in identifying features and properties of small ligands and their protein targets affecting the binding kinetics, but significant challenges include the long time scales involved in (un)binding events and the limited accuracy of empirical atomistic force fields (lacking, e.g., changes in electronic polarization). In an effort to overcome these hurdles, we propose a method that combines state-of-the-art enhanced sampling simulations and quantum mechanics/molecular mechanics (QM/MM) calculations at the BLYP/VDZ level to compute association free energy profiles and characterize the binding kinetics in terms of structure and dynamics of the transition state ensemble. We test our combined approach on the binding of the anticancer drug Imatinib to Src kinase, a well-characterized target for cancer therapy with a complex binding mechanism involving significant conformational changes. The results indicate significant changes in polarization along the binding pathways, which affect the predicted binding kinetics. This is likely to be of widespread importance in binding of ligands to protein targets.

Atomistic Analysis of ToxN and ToxI Complex Unbinding Mechanism

ToxIN is a triangular structure formed by three protein toxins (ToxNs) and three specific noncoding RNA antitoxins (ToxIs). To respond to stimuli, ToxI is preferentially degraded, releasing the ToxN. Thus, the dynamic character is essential in the normal function interactions between ToxN and ToxI. Here, equilibrated molecular dynamics (MD) simulations were performed to study the stability of ToxN and ToxI. The results indicate that ToxI adjusts the conformation of 3' and 5' termini to bind to ToxN. Steered molecular dynamics (SMD) simulations combined with the recently developed thermodynamic integration in 3nD (TI3nD) method were carried out to investigate ToxN unbinding from the ToxIN complex. The potentials of mean force (PMFs) and atomistic pictures suggest the unbinding mechanism as follows: (1) dissociation of the 5' terminus from ToxN, (2) missing the interactions involved in the 3' terminus of ToxI without three nucleotides (G31, A32, and A33), (3) starting to unfold for ToxI, (4) leaving the binding package of ToxN for three nucleotides of ToxI, (5) unfolding of ToxI. This work provides information on the structure-function relationship at the atomistic level, which is helpful for designing new potent antibacterial drugs in the future.

Exploring the PXR ligand binding mechanism with advanced Molecular Dynamics methods

The Pregnane X Receptor (PXR) is a ligand-activated transcription factor belonging to the nuclear receptor family. PXR can bind diverse drugs and environmental toxicants with different binding modes, making it an intriguing target for drug discovery. Here we investigated the binding mechanism of the SR12813 ligand to elucidate the significant steps, from the ligand entrance pathway into the binding cavity, to the ligand-induced conformational changes, and to the exploration of its alternative binding geometries. We used the advanced Molecular Dynamics-based methods implemented in the BiKi suite and developed specific methodological approaches to overcome the complexity induced by the buried and flexible binding cavity. The adopted methods provided a full dynamic description of the binding event and allowed rationalization of the observed multiple binding modes. These results suggest that the same approach could be exploited for the study of other binding processes with similar characteristics.

Machine learning accelerates MD-based binding pose prediction between ligands and proteins

Motivation

Fast and accurate prediction of protein–ligand binding structures is indispensable for structure-based drug design and accurate estimation of binding free energy of drug candidate molecules in drug discovery. Recently, accurate pose prediction methods based on short Molecular Dynamics (MD) simulations, such as MM-PBSA and MM-GBSA, among generated docking poses have been used. Since molecular structures obtained from MD simulation depend on the initial condition, taking the average over different initial conditions leads to better accuracy. Prediction accuracy of protein–ligand binding poses can be improved with multiple runs at different initial velocity.

Results

This paper shows that a machine learning method, called Best Arm Identification, can optimally control the number of MD runs for each binding pose. It allows us to identify a correct binding pose with a minimum number of total runs. Our experiment using three proteins and eight inhibitors showed that the computational cost can be reduced substantially without sacrificing accuracy. This method can be applied for controlling all kinds of molecular simulations to obtain best results under restricted computational resources.

Availability and implementation

Code and data are available on GitHub at https://github.com/tsudalab/bpbi.

Supplementary information

Supplementary data are available at Bioinformatics online.

Role of Resultant Dipole Moment in Mechanical Dissociation of Biological Complexes

Protein-peptide interactions play essential roles in many cellular processes and their structural characterization is the major focus of current experimental and theoretical research. Two decades ago, it was proposed to employ the steered molecular dynamics (SMD) to assess the strength of protein-peptide interactions. The idea behind using SMD simulations is that the mechanical stability can be used as a promising and an efficient alternative to computationally highly demanding estimation of binding affinity. However, mechanical stability defined as a peak in force-extension profile depends on the choice of the pulling direction. Here we propose an uncommon choice of the pulling direction along resultant dipole moment (RDM) vector, which has not been explored in SMD simulations so far. Using explicit solvent all-atom MD simulations, we apply SMD technique to probe mechanical resistance of ligand-receptor system pulled along two different vectors. A novel pulling direction—when ligand unbinds along the RDM vector—results in stronger forces compared to commonly used ligand unbinding along center of masses vector. Our observation that RDM is one of the factors influencing the mechanical stability of protein-peptide complex can be used to improve the ranking of binding affinities by using mechanical stability as an effective scoring function.

Binding Affinity via Docking: Fact and Fiction

In 1982, Kuntz et al. published an article with the title “A Geometric Approach to Macromolecule-Ligand Interactions”, where they described a method “to explore geometrically feasible alignment of ligands and receptors of known structure”. Since then, small molecule docking has been employed as a fast way to estimate the binding pose of a given compound within a specific target protein and also to predict binding affinity. Remarkably, the first docking method suggested by Kuntz and colleagues aimed to predict binding poses but very little was specified about binding affinity. This raises the question as to whether docking is the right tool to estimate binding affinity. The short answer is no, and this has been concluded in several comprehensive analyses. However, in this opinion paper we discuss several critical aspects that need to be reconsidered before a reliable binding affinity prediction through docking is realistic. These are not the only issues that need to be considered, but they are perhaps the most critical ones. We also consider that in spite of the huge efforts to enhance scoring functions, the accuracy of binding affinity predictions is perhaps only as good as it was 10–20 years ago. There are several underlying reasons for this poor performance and these are analyzed. In particular, we focus on the role of the solvent (water), the poor description of H-bonding and the lack of the systems’ true dynamics. We hope to provide readers with potential insights and tools to overcome the challenging issues related to binding affinity prediction via docking.

Simulated and experimental force spectroscopy of lysozyme on silica

The force spectra of proteins detaching from oxide surfaces measured by atomic force microscopy (AFM) often present complex patterns of peaks, which are difficult to correlate with individual bond-breaking events at the atomic scale. In this work we rationalize experimental AFM force spectra of hen-egg-white lysozyme detaching from silica by means of all-atom steered molecular dynamics (SMD) simulations. In particular, we demonstrate that the native tertiary structure of lysozyme is preserved if, and only if, its four intramolecular disulfide bridges are intact. Otherwise, the protein pulled off the surface undergoes severe unfolding, which is well captured by SMD simulations in explicit solvent. Implicit solvent simulations, on the contrary, wrongly predict protein unfolding even in the presence of S-S bridges, due to the lack of additional structural stabilization provided by the water's hydrogen-bond network within and surrounding the protein. On the basis of our combined experimental and theoretical findings, we infer that the rugged force spectra characteristic of lysozyme/silica interfaces are not due to the successive breaking of internal disulfide bonds leading to partial unfolding events. Rather, they reflect the detachment of several molecules bound to the same AFM tip, each anchored to the surface via multiple hydrogen and ionic bonds.

Computational modeling approaches to quantitative structure-binding kinetics relationships in drug discovery

Simple comparative correlation analyses and quantitative structure-kinetics relationship (QSKR) models highlight the interplay of kinetic rates and binding affinity as an essential feature in drug design and discovery. The choice of the molecular series, and their structural variations, used in QSKR modeling is fundamental to understanding the mechanistic implications of ligand and/or drug-target binding and/or unbinding processes. Here, we discuss the implications of linear correlations between kinetic rates and binding affinity constants and the relevance of the computational approaches to QSKR modeling.

Determining the Relative Binding Affinity of Ricin Toxin A Inhibitors by Using Molecular Docking and Nonequilibrium Work

Ricin is a ribosome-inactivating protein (RIP type 2) consisting of two subunits, ricin toxin A (RTA) and ricin toxin B (RTB). Because of its cytotoxicity, ricin has worried world authorities for its potential use as a chemical weapon; therefore, its inhibition is of great biotechnological interest. RTA is the target for inhibitor synthesis, and pterin derivatives are promising candidates to inhibit it. In this study, we used a combination of the molecular docking approach and fast steered molecular dynamics (SMD) to assess the correlation between nonequilibrium work, ⟨ W⟩, and the IC₅₀ for six RTA inhibitors. The results showed that molecular docking is a powerful tool to predict good bioactive poses of RTA inhibitors, and ⟨ W⟩ presented a strong correlation with IC₅₀ ( R² = 0.961). Such a profile ranked the RTA inhibitors better than the molecular docking approach. Therefore, the combination of docking and fast SMD simulation was shown to be a promising tool to distinguish RTA-active inhibitors from inactive ones and could be used as postdocking filtering approach.