Top Leads for Swine Influenza A/H1N1 Virus Revealed by Steered Molecular Dynamics Approach

Application of temperature-dependent and steered molecular dynamics simulation to screen anti-dengue compounds against Marburg virus

Marburg virus infections are extremely fatal with a fatality range of 23% to 90%, therefore there is an urgent requirement to design and develop efficient therapeutic molecules. Here, a comprehensive temperature-dependent molecular dynamics (MD) simulation method was implemented to identify the potential molecule from the anti-dengue compound library that can inhibit the function of the VP24 protein of Marburg. Virtual high throughput screening identified five effective binders of VP24 after screening 484 anti-dengue compounds. These compounds were treated in MD simulation at four different temperatures: 300, 340, 380, and 420 K. Higher temperatures showed dissociation of hit compounds from the protein. Further, triplicates of 100 ns MD simulation were conducted which showed that compounds ID = 118717693, and ID = 5361 showed strong stability with the protein molecule. These compounds were further validated using Δ G binding free energies and they showed: -30.38 kcal/mol, and -67.83 kcal/mol binding free energies, respectively. Later, these two compounds were used in steered MD simulation to detect its dissociation. Compound ID = 5361 showed the maximum pulling force of 199.02 kcal/mol/nm to dissociate the protein-ligand complex while ID = 118717693 had a pulling force of 101.11 kcal/mol/nm, respectively. This ligand highest number of hydrogen bonds with varying occupancies at 89.93%, 69.80%, 57.93%, 52.33%, and 50.63%. This study showed that ID = 5361 can bind with the VP24 strongly and has the potential to inhibit its function which can be validated in the in-vitro experiment.Communicated by Ramaswamy H. Sarma.

Calculation of binding affinity of JAK1 inhibitors via accurately computational estimation

Janus kinase 1 (JAK1) is a tyrosine kinase that is involved in the initiation of responses to a number of different cytokine receptor families. The JAK1-dependent pathway is a therapeutic target, and several JAK inhibitors have been developed thanks to intensive research. However, since the ATP binding sites of JAK family members are quite alike, JAK1 inhibitors can thus be less selective, resulting in unanticipated adverse effects. Despite this, minor variations in the ATP-binding site have been extensively used to find a variety of small compounds with different inhibitory properties. Stronger binding affinity of JAK1 inhibitors is believed to be able to reduce the negative effects, leading to better treatment results. Therefore, a thorough computational search that can effectively identify ligands with extremely high binding affinity for JAK1 to serve as promising inhibitors is required. Here, a method combining steered-molecular dynamic (SMD) simulations with a modified linear interaction energy (LIE) model has been developed to evaluate the binding affinities of known JAK1 inhibitors. The correlation coefficient between the estimated and experimental values was 0.72 and a root-mean-square error was 0.97 kcal•mol-1, revealing that the SMD/LIE method can precisely and quickly predict the binding free energies of JAK1 inhibitors. Furthermore, three marine fungus-derived compounds, namely hansforesters E, hansforesters G and tetroazolemycins B, were identified to be particularly promising JAK1 inhibitors, accordingly. These findings show that the SMD/LIE method has a lot of promise for in silico screening of possible JAK1 inhibitors from a vast number of compounds that are now accessible.Communicated by Ramaswamy H. Sarma.

Identifying inhibitors of NSP16-NSP10 of SARS-CoV-2 from large databases

The COVID-19 pandemic, which has already claimed millions of lives, continues to pose a serious threat to human health, requiring the development of new effective drugs. Non-structural proteins of SARS-CoV-2 play an important role in viral replication and infection. Among them, NSP16 (non-structured protein 16) and its cofactor NSP10 (non-structured protein 10) perform C2'-O methylation at the 5' end of the viral RNA, which promotes efficient virus replication. Therefore, the NSP16-NSP10 complex becomes an attractive target for drug development. Using a multi-step virtual screening protocol which includes Lipinski's rule, docking, steered molecular dynamics and umbrella sampling, we searched for potential inhibitors from the PubChem and anti-HIV databases. It has been shown that CID 135566620 compound from PubChem is the best candidate with an inhibition constant in the sub-μM range. The Van der Waals interaction was found to be more important than the electrostatic interaction in the binding affinity of this compound to NSP16-NSP10. Further in vitro and in vivo studies are needed to test the activity of the identified compound against COVID-19.Communicated by Ramaswamy H. Sarma.

Steered Molecular Dynamics Simulations Study on FABP4 Inhibitors

Ordinary small molecule de novo drug design is time-consuming and expensive. Recently, computational tools were employed and proved their efficacy in accelerating the overall drug design process. Molecular dynamics (MD) simulations and a derivative of MD, steered molecular dynamics (SMD), turned out to be promising rational drug design tools. In this paper, we report the first application of SMD to evaluate the binding properties of small molecules toward FABP4, considering our recent interest in inhibiting fatty acid binding protein 4 (FABP4). FABP4 inhibitors (FABP4is) are small molecules of therapeutic interest, and ongoing clinical studies indicate that they are promising for treating cancer and other diseases such as metabolic syndrome and diabetes.

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.

Antibody-nanobody combination increases their neutralizing activity against SARS-CoV-2 and nanobody H11-H4 is effective against Alpha, Kappa and Delta variants

The global spread of COVID-19 is devastating health systems and economies worldwide. While the use of vaccines has yielded encouraging results, the emergence of new variants of SARS-CoV-2 shows that combating COVID-19 remains a big challenge. One of the most promising treatments is the use of not only antibodies, but also nanobodies. Recent experimental studies revealed that the combination of antibody and nanobody can significantly improve their neutralizing ability through binding to the SARS-CoV-2 spike protein, but the molecular mechanisms underlying this observation remain largely unknown. In this work, we investigated the binding affinity of the CR3022 antibody and H11-H4 nanobody to the SARS-CoV-2 receptor binding domain (RBD) using molecular modeling. Both all-atom steered molecular dynamics simulations and coarse-grained umbrella sampling showed that, consistent with the experiment, CR3022 associates with RBD more strongly than H11-H4. We predict that the combination of CR3022 and H11-H4 considerably increases their binding affinity to the spike protein. The electrostatic interaction was found to control the association strength of CR3022, but the van der Waals interaction dominates in the case of H11-H4. However, our study for a larger set of nanobodies and antibodies showed that the relative role of these interactions depends on the specific complex. Importantly, we showed Beta, Gamma, Lambda, and Mu variants reduce the H11-H4 activity while Alpha, Kappa and Delta variants increase its neutralizing ability, which is in line with experiment reporting that the nanobody elicited from the llama is very promising for fighting against the Delta variant.

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.

Structural analysis of TrkA mutations in patients with congenital insensitivity to pain reveals PLCγ as an analgesic drug target

Chronic pain is a major health issue, and the search for new analgesics has become increasingly important because of the addictive properties and unwanted side effects of opioids. To explore potentially new drug targets, we investigated mutations in the NTRK1 gene found in individuals with congenital insensitivity to pain with anhidrosis (CIPA). NTRK1 encodes tropomyosin receptor kinase A (TrkA), the receptor for nerve growth factor (NGF) and that contributes to nociception. Molecular modeling and biochemical analysis identified mutations that decreased the interaction between TrkA and one of its substrates and signaling effectors, phospholipase Cγ (PLCγ). We developed a cell-permeable phosphopeptide derived from TrkA (TAT-pQYP) that bound the Src homology domain 2 (SH2) of PLCγ. In HEK-293T cells, TAT-pQYP inhibited the binding of heterologously expressed TrkA to PLCγ and decreased NGF-induced, TrkA-mediated PLCγ activation and signaling. In mice, intraplantar administration of TAT-pQYP decreased mechanical sensitivity in an inflammatory pain model, suggesting that targeting this interaction may be analgesic. The findings demonstrate a strategy to identify new targets for pain relief by analyzing the signaling pathways that are perturbed in CIPA.

Umbrella Sampling-Based Method to Compute Ligand-Binding Affinity

Many proteins have a solvent-exposed binding cleft, which permits their inhibitors to bind and unbind without significant protein conformation transforms. The binding/unbinding pathways of these protein-inhibitor complexes can be rather straightforwardly sampled by using umbrella sampling (US) simulation methods. During a US simulation, the C_α atoms of the protein are restrained via a harmonic force. The potential of mean force (PMF) along the binding pathway can be estimated by using the weighted histogram analysis method (WHAM). The binding affinity is then computed as the difference in PMF between the binding and unbinding states.

The effects of mutation on the drug binding affinity of Neuraminidase: case study of Capsaicin using steered molecular dynamics simulation

The influenza virus is an important respiratory pathogen that causes many incidences of diseases and even death each year. One of the primary factors of this virus is the Neuraminidase surface protein, which causes the virus to leave the host cell and spread to new target cells. The main antiviral medication for influenza is designed as a protein inhibitor ligand that prevents further spread of the disease, and eventually relieves the emerged symptoms. The effectiveness of such inhibitory drugs is highly associated with their binding affinity. In this paper, the binding affinity of an herbal ligand of Capsaicin bound to Neuraminidase of the influenza virus is investigated using steered molecular dynamics (SMD) simulation. Since mutations of the virus directly impact the binding affinity of the inhibitory drugs, different mutations were generated by using Mutagenesis module. The rapid spread of infection during the avian influenza A/H5N1 epidemic has raised concerns about far more dangerous consequences if the virus becomes resistant to current drugs. Currently, oseltamivir (Tamiflu), zanamivir (Relenza), pramivir (Rapivab), and laninamivir (Inavir) are increasingly used to treat the flu. However, with the rapid evolution of the virus, some drug-resistant strains are emerging. Therefore, it is very important to seek alternative therapies and identify the roots of drug resistance. Obtained results demonstrated a reduced binding affinity for the applied mutations. This reduction in binding affinity will cause the virus mutation to become resistant to the drug, which will spread the disease and make it more difficult to treat. From a molecular prospect, this decrease in binding affinity is due to the loss of a number of effective bonds between the ligand and the receptor, which occurs with mutations of the wild-type (WT) species. The results of the present study can be used in the rational design of novel drugs that are compatible with specific mutations.

Searching and designing potential inhibitors for SARS-CoV-2 Mpro from natural sources using atomistic and deep-learning calculations

The spread of severe acute respiratory syndrome coronavirus 2 novel coronavirus (SARS-CoV-2) worldwide has caused the coronavirus disease 2019 (COVID-19) pandemic. A hundred million people were infected, resulting in several millions of death worldwide. In order to prevent viral replication, scientists have been aiming to prevent the biological activity of the SARS-CoV-2 main protease (3CL pro or Mpro). In this work, we demonstrate that using a reasonable combination of deep-learning calculations and atomistic simulations could lead to a new approach for developing SARS-CoV-2 main protease (Mpro) inhibitors. Initially, the binding affinities of the natural compounds to SARS-CoV-2 Mpro were estimated via atomistic simulations. The compound tomatine, thevetine, and tribuloside could bind to SARS-CoV-2 Mpro with nanomolar/high-nanomolar affinities. Secondly, the deep-learning (DL) calculations were performed to chemically alter the top-lead natural compounds to improve ligand-binding affinity. The obtained results were then validated by free energy calculations using atomistic simulations. The outcome of the research will probably boost COVID-19 therapy.

The hybrid DeepFrag/atomistic simulation approach could lead to a new scheme for developing SARS-CoV-2 3CLpro/Mpro inhibitors.

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.

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.

Binding of inhibitors to the monomeric and dimeric SARS-CoV-2 Mpro

SARS-CoV-2 rapidly infects millions of people worldwide since December 2019. There is still no effective treatment for the virus, resulting in the death of more than one million patients. Inhibiting the activity of SARS-CoV-2 main protease (Mpro), 3C-like protease (3CLP), is able to block the viral replication and proliferation. In this context, our study has revealed that in silico screening for inhibitors of SARS-CoV-2 Mpro can be reliably done using the monomeric structure of the Mpro instead of the dimeric one. Docking and fast pulling of ligand (FPL) simulations for both monomeric and dimeric forms correlate well with the corresponding experimental binding affinity data of 24 compounds. The obtained results were also confirmed via binding pose and noncovalent contact analyses. Our study results show that it is possible to speed up computer-aided drug design for SARS-CoV-2 Mpro by focusing on the monomeric form instead of the larger dimeric one.

Computational Determination of Potential Inhibitors of SARS-CoV-2 Main Protease

The novel coronavirus (SARS-CoV-2) has infected several million people and caused thousands of deaths worldwide since December 2019. As the disease is spreading rapidly all over the world, it is urgent to find effective drugs to treat the virus. The main protease (Mpro) of SARS-CoV-2 is one of the potential drug targets. Therefore, in this context, we used rigorous computational methods, including molecular docking, fast pulling of ligand (FPL), and free energy perturbation (FEP), to investigate potential inhibitors of SARS-CoV-2 Mpro. We first tested our approach with three reported inhibitors of SARS-CoV-2 Mpro, and our computational results are in good agreement with the respective experimental data. Subsequently, we applied our approach on a database of ∼4600 natural compounds, as well as 8 available HIV-1 protease (PR) inhibitors and an aza-peptide epoxide. Molecular docking resulted in a short list of 35 natural compounds, which was subsequently refined using the FPL scheme. FPL simulations resulted in five potential inhibitors, including three natural compounds and two available HIV-1 PR inhibitors. Finally, FEP, the most accurate and precise method, was used to determine the absolute binding free energy of these five compounds. FEP results indicate that two natural compounds, cannabisin A and isoacteoside, and an HIV-1 PR inhibitor, darunavir, exhibit a large binding free energy to SARS-CoV-2 Mpro, which is larger than that of 13b, the most reliable SARS-CoV-2 Mpro inhibitor recently reported. The binding free energy largely arises from van der Waals interaction. We also found that Glu166 forms H-bonds to all of the inhibitors. Replacing Glu166 by an alanine residue leads to ∼2.0 kcal/mol decreases in the affinity of darunavir to SARS-CoV-2 Mpro. Our results could contribute to the development of potential drugs inhibiting SARS-CoV-2.

Investigation on penetration of saffron components through lipid bilayer bound to spike protein of SARS-CoV-2 using steered molecular dynamics simulation

A coronavirus identified as COVID-19 is the reason for an infection outbreak which is started in December 2019. NO completely effective drugs and treatments are not recognized for this virus. Recently, saffron and its compounds were used to treat different viral diseases. Saffron extract and its major ingredients have shown antiviral effects. In this study, the steered molecular dynamics simulation was used for investigating the effect of four main components of saffron that include: crocin, crocetin, safranal, and picrocrocin as candidate for drug molecules, on COVID-19. The binding energies between drug molecules and spike protein and the main protease of the virus were evaluated. The obtained results based on Lennard-Jones and electrostatic potentials demonstrated that crocetin has a high affinity towards spike protein and also the main protease of the virus. Also, the quantum mechanics calculations elucidated that the crocetin could overcome energy barrier of lipid bilayer with strong dipole moment and polarizability. The pharmacokinetic and ADMET properties proved that crocetin could be a suitable drug candidate. So, crocetin could be a promising drug for treatment of COVID-19.

Estimation of the ligand-binding free energy of checkpoint kinase 1 via non-equilibrium MD simulations

Checkpoint kinase 1 (CHK1) is a serine/threonine-protein kinase that is involved in cell cycle regulation in eukaryotes. Inhibition of CHK1 is thus considered as a promising approach in cancer therapy. In this study, the fast pulling of ligand (FPL) process was applied to predict the relative binding affinities of CHK1 inhibitors using non-equilibrium molecular dynamics (MD) simulations. The work of external harmonic forces to pull the ligand out of the binding cavity strongly correlated with the experimental binding affinity of CHK1 inhibitors with the correlation coefficient of R = -0.88 and an overall root mean square error (RMSE) of 0.99 kcal/mol. The data indicate that the FPL method is highly accurate in predicting the relative binding free energies of CHK1 inhibitors with an affordable CPU time. A new set of molecules were designed based on the molecular modeling of interactions between the known inhibitor and CHK1 as inhibitory candidates. Molecular docking and FPL results exhibited that the binding affinities of developed ligands were similar to the known inhibitor in interaction with the catalytic site of CHK1, producing very potential CHK1 inhibitors of that the inhibitory activities should be further evaluated in vitro.

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.

Fast Estimation of the Blood-Brain Barrier Permeability by Pulling a Ligand through a Lipid Membrane

The blood–brain barrier (BBB) is a physical barrier that regulates the homeostasis of the neural microenvironment. A relative estimate of the BBB permeability, which is important for drug design, may be experimentally provided by the logBB (the blood–brain concentration ratio) and the logPS (permeability–surface-area product), while many computational methods aim to identify key properties that correlate well with these quantities. Although currently existing computational methods (e.g., quantitative structure activity relation) have made a significant contribution in screening various compounds that could potentially translocate through the BBB, they are unable to provide a physical explanation of the underlying processes and they can often be computationally demanding. Here, we use steered molecular dynamics simulation to estimate the BBB permeability of various compounds on the basis of simple lipid–membrane models by computing the nonequilibrium work, W_neq, produced by pulling the compounds through the membrane. We found that the values of W_neq correlate remarkably well with logBB and logPS for a range of compounds and different membrane types and pulling speeds, independently of the choice of force field. Moreover, our results provide insight into the role of hydrogen bonds, the energetic barriers, and the forces exerted on the ligands during their pulling. Our method is computationally easy to implement and fast. Therefore, we anticipate that it could provide a reliable prescreening tool for estimating the relative permeability of the BBB to various substances.

Calculating the absolute binding free energy of the insulin dimer in an explicit solvent

Insulin is a significant hormone in the regulation of glucose level in the blood. Its monomers bind to each other to form dimers or hexamers through a complex process. To study the binding of the insulin dimer, we first calculate its absolute binding free energy by the steered molecular dynamics method and the confinement method based on a fictitious thermodynamic cycle. After considering some special correction terms, the final calculated binding free energy at 298 K is −8.97 ± 1.41 kcal mol⁻¹, which is close to the experimental value of −7.2 ± 0.8 kcal mol⁻¹. Furthermore, we discuss the important residue–residue interactions between the insulin monomers, including hydrophobic interactions, π–π interactions and hydrogen bond interactions. The analysis reveals five key residues, Vla^B12, Tyr^B16, Phe^B24, Phe^B25, and Tyr^B26, for the dimerization of the insulin. We also perform MM-PBSA calculations for the wild-type dimer and some mutants and study the roles of the key residues by the change of the binding energy of the insulin dimer.

In this paper, we calculate the absolute binding free energy of an insulin dimer by steered MD method. The result of −8.97 kcal mol⁻¹ is close to the experimental value −7.2 kcal mol⁻¹. We also analyze the residue–residue interactions.

Structural Analysis of PfSec62-Autophagy Interacting Motifs (AIM) and PfAtg8 Interactions for Its Implications in RecovER-phagy in Plasmodium falciparum

Autophagy is a degradative pathway associated with many pathological and physiological processes crucial for cell survival. During ER stress, while selective autophagy occurs via ER-phagy, the re-establishment of physiologic ER homeostasis upon resolution of a transient ER stress is mediated by recovER-phagy. Recent studies demonstrated that recovER-phagy is governed via association of Sec62 as an ER-resident autophagy receptor through its autophagy interacting motifs (AIM)/LC3-interacting region (LIR) toAtg8/LC3. Atg8 is an autophagy protein, which is central to autophagosome formation and maturation. Plasmodium falciparum Atg8 (PfAtg8) has both autophagic and non-autophagic functions critical for parasite survival. Since Plasmodium also has Sec62 in the ER membrane and is prone to ER stress due to drastic transformation during their complex intraerythrocytic cycle; hence, we initiated the studies to check whether recovER-phagy occurs in the parasite. To achieve this, a comprehensive study based on the computational approaches was carried out. This study embarks upon identification of AIM sequences in PfSec62 by carrying out peptide-protein docking simulations and comparing the interactions of these AIMs with PfAtg8, based on the molecular dynamic simulations. Detailed analysis is based on electrostatic surface complementarity, peptide-protein interaction strength, mapping of non-covalent bond interactions and rupture force calculated from steered MD simulations. Potential mean forces and unbinding free energies (ΔG_dissociation) using Jarzynski's equality were also computed for the AIM/LIR motif complexes with PfAtg8/HsLC3 autophagy proteins to understand their dissociation free energy profiles and thereby their binding affinities and stability of the peptide-protein complexes. Through this study, we predict Sec62 mediated recovER-phagy in Plasmodium falciparum, which might open new avenues to explore novel drug targets for antimalarial drug discovery.

Effective Estimation of Ligand-Binding Affinity Using Biased Sampling Method

The binding between two biomolecules is one of the most critical factors controlling many bioprocesses. Therefore, it is of great interest to derive a reliable method to calculate the free binding energy between two biomolecules. In this work, we have demonstrated that the binding affinity of ligands to proteins can be determined through biased sampling simulations. The umbrella sampling (US) method was applied on 20 protein-ligand complexes, including the cathepsin K (CTSK), type II dehydroquinase (DHQase), heat shock protein 90 (HSP90), and factor Xa (FXa) systems. The ligand-binding affinity was evaluated as the difference between the largest and smallest values of the free-energy curve, which was obtained via a potential of mean force analysis. The calculated affinities differ sizably from the previously reported experimental values, with an average difference of ∼3.14 kcal/mol. However, the calculated results are in good correlation with the experimental data, with correlation coefficients of 0.76, 0.87, 0.96, and 0.97 for CTSK, DHQase, HSP90, and FXa, respectively. Thus, the binding free energy of a new ligand can be reliably estimated using our US approach. Furthermore, the root-mean-square errors (RMSEs) of binding affinity of these systems are 1.13, 0.90, 0.37, and 0.25 kcal/mol, for CTSK, DHQase, HSP90, and FXa, respectively. The small RMSE values indicate the good precision of the biased sampling method that can distinguish the ligands exhibiting similar binding affinities.

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.

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.

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.

Probing the Binding Affinity by Jarzynski's Nonequilibrium Binding Free Energy and Rupture Time

Binding affinity of a small ligand to a receptor is the important quantity in drug design, and it might be characterized by different quantities. The most popular one is the binding free energy, which can be estimated by several methods in conventional molecular dynamics simulation. So far in steered molecular dynamics (SMD), one can use either the rupture force or nonequilibrium pulling work as a measure for binding affinity. In this paper, we have shown that the nonequilibrium binding free energy Δ G_neq^Jar, obtained by Jarzynski's equality at a finite pulling speed, has good correlation with experimental data on inhibition constants, implying that this quantity can be used as a good scoring function for binding affinity. A similar correlation has also been disclosed for binding and unbinding free energy barriers. Applying the SMD method to unbinding of 23 small compounds from the binding site of β-lactamase protein, a bacteria-produced enzyme, we have demonstrated that the rupture or unbinding time strongly correlates with experimental data with correlation level R ≈ 0.84. As follows from the Jarzynski's equality, the rupture time depends on the unbinding barrier exponentially. We show that Δ G_neq^Jar, the rupture time, and binding and unbinding free energy barriers are good descriptors for binding affinity. Our observation may be useful for fast screening of potential leads as the SMD simulation is not time-consuming. On the basis of nonequilibrium simulation, we disclosed that, in agreement with the experiment, the binding time is much longer than the unbinding one.

Single-molecule techniques in biophysics: a review of the progress in methods and applications

Single-molecule biophysics has transformed our understanding of biology, but also of the physics of life. More exotic than simple soft matter, biomatter lives far from thermal equilibrium, covering multiple lengths from the nanoscale of single molecules to up to several orders of magnitude higher in cells, tissues and organisms. Biomolecules are often characterized by underlying instability: multiple metastable free energy states exist, separated by levels of just a few multiples of the thermal energy scale k _B T, where k _B is the Boltzmann constant and T absolute temperature, implying complex inter-conversion kinetics in the relatively hot, wet environment of active biological matter. A key benefit of single-molecule biophysics techniques is their ability to probe heterogeneity of free energy states across a molecular population, too challenging in general for conventional ensemble average approaches. Parallel developments in experimental and computational techniques have catalysed the birth of multiplexed, correlative techniques to tackle previously intractable biological questions. Experimentally, progress has been driven by improvements in sensitivity and speed of detectors, and the stability and efficiency of light sources, probes and microfluidics. We discuss the motivation and requirements for these recent experiments, including the underpinning mathematics. These methods are broadly divided into tools which detect molecules and those which manipulate them. For the former we discuss the progress of super-resolution microscopy, transformative for addressing many longstanding questions in the life sciences, and for the latter we include progress in 'force spectroscopy' techniques that mechanically perturb molecules. We also consider in silico progress of single-molecule computational physics, and how simulation and experimentation may be drawn together to give a more complete understanding. Increasingly, combinatorial techniques are now used, including correlative atomic force microscopy and fluorescence imaging, to probe questions closer to native physiological behaviour. We identify the trade-offs, limitations and applications of these techniques, and discuss exciting new directions.

Evaluation of the absolute affinity of neuraminidase inhibitor using steered molecular dynamics simulations

The absolute free energy difference of binding (ΔG) between neuraminidase and its inhibitor was evaluated using fast pulling of ligand (FPL) method over steered molecular dynamics (SMD) simulations. The metric was computed through linear interaction approximation. Binding nature was described by free energy differences of electrostatic and van der Waals (vdW) interactions. The finding indicates that vdW metric is dominant over electrostatics in binding process. The computed values are in good agreement with experimental data with a correlation coefficient of R=0.82 and error of σΔG_exp=2.2kcal/mol. The results were observed using Amber99SB-ILDN force field in comparison with CHARMM27 and GROMOS96 43a1 force fields. Obtained results may stimulate the search for an Influenza therapy.

Protocol for fast screening of multi-target drug candidates: Application to Alzheimer's disease

The treatment of many diseases may require drugs that are capable to attack multiple targets simultaneously. Obviously, the virtual screening of multi-target drug candidates is much more time consuming compared to the single-target case. This, in particular, concerns the last step of virtual screening where the binding free energy is computed by conventional molecular dynamics simulation. To overcome this difficulty we propose a simple protocol which is relied on the fast steered molecular dynamics simulation and on available experimental data on binding affinity of reference ligand to a given target. Namely, first we compute non-equilibrium works generated during pulling ligands from the binding site using the steered molecular dynamics method. Then as top leads we choose only those compounds that have the non-equilibrium work larger than that of a reference compound for which the binding free energy has been already known from experiment. Despite many efforts no cures for AD (Alzheimer's disease) have been found. One of possible reasons for this failure is that drug candidates were developed for a single target, while there are exist many possible pathways to AD. Applying our new protocol to five targets including amyloid beta fibril, peroxisome proliferator-activated receptor γ, retinoic X receptor α, β- and γ-secretases, we have found two potential drugs (CID 16040294 and CID 9998128) for AD from the large PubChem database. We have also shown that these two ligands can interfere with the activity of popular Acetylcholinesterase target through strong binding towards it.

Evaluation of protein-ligand affinity prediction using steered molecular dynamics simulations

In computational drug design, ranking a series of compound analogs in a manner that is consistent with experimental affinities remains a challenge. In this study, we evaluated the prediction of protein-ligand binding affinities using steered molecular dynamics simulations. First, we investigated the appropriate conditions for accurate predictions in these simulations. A conic harmonic restraint was applied to the system for efficient sampling of work values on the ligand unbinding pathway. We found that pulling velocity significantly influenced affinity predictions, but that the number of collectable trajectories was less influential. We identified the appropriate pulling velocity and collectable trajectories for binding affinity predictions as 1.25 Å/ns and 100, respectively, and these parameters were used to evaluate three target proteins (FK506 binding protein, trypsin, and cyclin-dependent kinase 2). For these proteins using our parameters, the accuracy of affinity prediction was higher and more stable when Jarzynski's equality was employed compared with the second-order cumulant expansion equation of Jarzynski's equality. Our results showed that steered molecular dynamics simulations are effective for predicting the rank order of ligands; thus, they are a potential tool for compound selection in hit-to-lead and lead optimization processes.

Discovery of DNA dyes Hoechst 34580 and 33342 as good candidates for inhibiting amyloid beta formation: in silico and in vitro study

Combining Lipinski's rule with the docking and steered molecular dynamics simulations and using the PubChem data base of about 1.4 million compounds, we have obtained DNA dyes Hoechst 34580 and Hoechst 33342 as top-leads for the Alzheimer's disease. The binding properties of these ligands to amyloid beta (Aβ) fibril were thoroughly studied by in silico and in vitro experiments. Hoechst 34580 and Hoechst 33342 prefer to locate near hydrophobic regions with binding affinity mainly governed by the van der Waals interaction. By the Thioflavin T assay, it was found that the inhibition constant IC50 ≈ 0.86 and 0.68 μM for Hoechst 34580 and Hoechst 33342, respectively. This result qualitatively agrees with the binding free energy estimated using the molecular mechanic-Poisson Boltzmann surface area method and all-atom simulations with the AMBER-f99SB-ILDN force field and water model TIP3P. In addition, DNA dyes have the high capability to cross the blood brain barrier. Thus, both in silico and in vitro experiments have shown that Hoechst 34580 and 33342 are good candidates for treating the Alzheimer's disease by inhibiting Aβ formation.

Screening of hepatitis C NS5B polymerase inhibitors containing benzothiadiazine core: a steered molecular dynamics

Hepatic C virus (HCV) is a global health problem, resulting in liver cirrhosis and inflammation that can develop to hepatocellular carcinoma and fatality. The NS5B polymerase of HCV plays an important role in viral RNA replication process, making it an attractive therapeutic target for design and development of anti-HCV drugs. To search new potent compounds against the HCV NS5B polymerase, the molecular docking and the steered molecular dynamics (SMD) simulation techniques were performed. The potential potent inhibitors of the NS5B polymerase were screened out from the ZINC database using structural similarity search and molecular docking technique. Five top-hit compounds (the ZINC compounds 49888724, 49054741, 49777239, 49793673, and 49780355) were then studied by the SMD simulations based on the hypothesis that a high rupture force relates to a high binding efficiency. The results demonstrated that the ZINC compound 49888724 had a greater maximum rupture force, reflecting a good binding strength and inhibitory potency than known inhibitors and the rest four ZINC compounds. Therefore, our finding indicated that the ZINC compound 49888724 is a potential candidate to be a novel NS5B inhibitor for further design. Besides, the van der Waals interaction could be considered as the main contribution for stabilizing the NS5B-ligand complex.

Ligand binding to anti-cancer target CD44 investigated by molecular simulations

CD44 is a cell-surface glycoprotein and receptor for hyaluronan, one of the major components of the tumor extracellular matrix. There is evidence that the interaction between CD44 and hyaluronan promotes breast cancer metastasis. Recently, the molecule F-19848A was shown to inhibit hyaluronan binding to receptor CD44 in a cell-based assay. In this study, we investigated the mechanism and energetics of F-19848A binding to CD44 using molecular simulation. Using the molecular mechanics/Poisson Boltzmann surface area (MM-PBSA) method, we obtained the binding free energy and inhibition constant of the complex. The van der Waals (vdW) interaction and the extended portion of F-19848A play key roles in the binding affinity. We screened natural products from a traditional Chinese medicine database to search for CD44 inhibitors. From combining pharmaceutical requirements with docking and molecular dynamics simulations, we found ten compounds that are potentially better or equal to the F-19848A ligand at binding to CD44 receptor. Therefore, we have identified new candidates of CD44 inhibitors, based on molecular simulation, which may be effective small molecules for the therapy of breast cancer.

A New Method for Navigating Optimal Direction for Pulling Ligand from Binding Pocket: Application to Ranking Binding Affinity by Steered Molecular Dynamics

In this paper we present a new method for finding the optimal path for pulling a ligand from the binding pocket using steered molecular dynamics (SMD). Scoring function is defined as the steric hindrance caused by a receptor to ligand movement. Then the optimal path corresponds to the minimum of this scoring function. We call the new method MSH (Minimal Steric Hindrance). Contrary to existing navigation methods, our approach takes into account the geometry of the ligand while other methods including CAVER only consider the ligand as a sphere with a given radius. Using three different target + receptor sets, we have shown that the rupture force Fmax and nonequilibrium work Wpull obtained based on the MSH method show a much higher correlation with experimental data on binding free energies compared to CAVER. Furthermore, Wpull was found to be a better indicator for binding affinity than Fmax. Thus, the new MSH method is a reliable tool for obtaining the best direction for ligand exiting from the binding site. Its combination with the standard SMD technique can provide reasonable results for ranking binding affinities using Wpull as a scoring function.