Rapid and Accurate Prediction of Binding Free Energies for Saquinavir-Bound HIV-1 Proteases

A Study on the Binding Mechanism and the Impact of Key Residue Mutations between SND1 and MTDH Peptide through Molecular Dynamics Simulations

Metastasis of breast cancer is the main cause of death for patients with breast cancer. The interaction between metadherin (MTDH) and staphylococcal nuclease domain 1 (SND1) plays a pivotal role in promoting breast cancer development. However, the binding details between MTDH and SND1 remain unclear. In this study, we employed all-atom molecular dynamics simulations (MDs) and conducted binding energy calculations to investigate the binding details and the impact of key residue mutations on binding. The mutations in key residues have not significantly affected the overall stability of the structure and the fluctuation of residues near the binding site; they have exerted a substantial impact on the binding of SND1 and MTDH peptide. The electrostatic interactions and van der Waals interactions play an important role in the binding of SND1 and the MTDH peptide. The mutations in the key residues have a significant impact on electrostatic and van der Waals interactions, resulting in weakened binding. The energy contributions of key residues mainly come from the electrostatic energy and van der Waals interactions of the side chain. In addition, the key residues form an intricate and stable network of hydrogen bonds and salt-bridge interactions with the MTDH peptide. The mutations in key residues have directly disrupt the interactions formed between SND1 and MTDH peptide, consequently leading to changes in the binding mode of the MTDH peptide. These analyses unveil the detailed atomic-level interaction mechanism between SND1 and the MTDH peptide, providing a molecular foundation for the development of antibreast cancer drugs.

Revealing the drug resistance mechanism of saquinavir due to G48V and V82F mutations in subtype CRF01_AE HIV-1 protease: molecular dynamics simulation and binding free energy calculations

Human immunodeficiency virus-1 (HIV-1) protease is one of the important targets in AIDS therapy. The majority of HIV infections are caused due to non-B subtypes in developing countries. The co-occurrence of mutations along with naturally occurring polymorphisms in HIV-1 protease cause resistance to the FDA approved drugs, thereby posing a major challenge in the treatment of antiretroviral therapy. In this work, the resistance mechanism against SQV due to active site mutations G48V and V82F in CRF01_AE (AE) protease was explored. The binding free energy calculations showed that the direct substitution of valine at position 48 introduces a bulkier side chain, directly impairing the interaction with SQV in the binding pocket. Also, the intramolecular hydrogen bonding network of the neighboring residues is altered, indirectly affecting the binding of SQV. Interestingly, the substitution of phenylalanine at position 82 induces conformational changes in the 80's loop and the flap region, thereby favoring the binding of SQV. The V82F mutant structure also maintains similar intramolecular hydrogen bond interactions as observed in AE-WT.Communicated by Ramaswamy H. Sarma.

Functionalized carbon nanotubes as an alternative to traditional anti-HIV-1 protease inhibitors: An understanding towards Nano-medicine development through MD simulations

The Human Immunodeficiency Virus (HIV) has been the source of epidemic infection of AIDS for a longer period. One of the most difficult tasks is identifying novel medications that can help to decrease or control this global health hazard by overcoming drug resistance. In recent decades' nanoparticles are emerging as extremely relevant in drug delivery platforms. In the current study, the pristine (SWCNT) and hydroxyl functionalized (SWCNT-OH) versions of the SWCNT were investigated as inhibitors against the wild-type (WT) and three key mutants of HIV-1 protease (HIV-pr) (I50V, V82A, and I84V). Molecular docking of SWCNT in the catalytic domain and running all-atom MD simulations of all complexes are also part of this project. A thorough inspection of conformational dynamics from 50 ns trajectories reveals that both the pristine and SWCNT-OH can fit right to the pocket region of HIV-pr and govern flap dynamics. The binding affinity of the four HIV-pr-SWCNT/SWCNT-OH complexes was further investigated using MM-PBSA-dependent binding free energy studies. In most mutants and WT systems, SWCNT-OH was reported to bind proportionately many folds (kcal/mol) more than pristine SWCNTs. Hence, SWCNTs are possible HIV-pr inhibitors in terms of their stable existence in the pocket area, stronger binding to the protease, and regulation of flap dynamics in controlling the active site volume, which have vast potential for applications against drug resistance.

Classical and Machine Learning Methods for Protein - Ligand Binding Free Energy Estimation

Binding free energy estimation of drug candidates to their biomolecular target is one of the best quantitative estimators in computer-aided drug discovery. Accurate binding free energy estimation is still a challengeable task even after decades of research, along with the complexity of the algorithm, time-consuming procedures, and reproducibility issues. In this review, we have discussed the advantages and disadvantages of diverse free energy methods like Thermodynamic Integration (TI), Bennett's Acceptance Ratio (BAR), Free Energy Perturbation (FEP), alchemical methods. Moreover, we discussed the possible application of the machine learning method in protein-ligand binding free energy estimation.

Effect of Stapling on the Thermodynamics of mdm2-p53 Binding

Protein-protein interaction (PPI) is one of the key regulatory features driving biomolecular processes and hence is targeted for designing therapeutics against diseases. Small peptides are a new and emerging class of therapeutics owing to their high specificity and low toxicity. For achieving efficient targeting of the PPI, amino acid side chains are often stapled together, resulting in the rigidification of these peptides. Exploring the scope of these peptides demands a comprehensive understanding of their working principle. In this work, two stapled p53 peptides have been considered to delineate their binding mechanism with mdm2 using computational approaches. The addition of stapling agent protects the secondary structure of the peptides even in the case of thermal and chemical denaturation. Although the introduction of a stapling agent increases the hydrophobicity of the peptide, the enthalpic stabilization decreases. This is overcome by the lowering of the entropic penalty, and the overall binding affinity improves. The mechanistic insights into the benefit of peptide stapling can be adopted for further improvement of peptide therapeutics.

Computational Design of PDZ-Peptide Binding

This chapter describes two computational methods for PDZ-peptide binding: high-throughput computational protein design (CPD) and a medium-throughput approach combining molecular dynamics for conformational sampling with a Poisson-Boltzmann (PB) Linear Interaction Energy for scoring. A new CPD method is outlined, which uses adaptive Monte Carlo simulations to efficiently sample peptide variants that tightly bind a PDZ domain, and provides at the same time precise estimates of their relative binding free energies. A detailed protocol is described based on the Proteus CPD software. The medium-throughput approach can be performed with standard MD and PB software, such as NAMD and Charmm. For 40 complexes between Tiam1 and peptide ligands, it gave high a2ccuracy, with mean errors of around 0.5 kcal/mol for relative binding free energies and no large errors. It requires a moderate amount of parameter fitting before it can be applied, and its transferability to other protein families is still untested.

Ca2+-based allosteric switches and shape shifting in RGLG1 VWA domain

RGLG1 is an E3 ubiquitin ligase in Arabidopsis thaliana that participates in ABA signaling and regulates apical dominance. Here, we present crystal structures of RGLG1 VWA domain, revealing two novel calcium ions binding sites (NCBS1 and NCBS2). Furthermore, the structures with guided mutagenesis in NCBS1 prove that Ca²⁺ ions play important roles in controlling conformational change of VWA, which is stabilized in open state with Ca²⁺ bound and converted to closed state after Ca²⁺ removal. This allosteric regulation mechanism is distinct from the ever reported one involving the C-terminal helix in integrin α and β I domains. The mutation of a key residue in NCBS2 do not abolish its Ca²⁺-binding potential, with no conformational change. MD simulations reveals that open state of RGLG1 VWA has higher ligand affinity than its closed state, consisting with integrin. Structural comparison of ion-free-MIDAS with Mg²⁺-MIDAS reveals that Mg²⁺ binding to MIDAS does not induce conformational change. With acquisition of first structure of plant VWA domain in both open state and closed state, we carefully analyze the conformational change and propose a totally new paradigm for its transition of open-closed states, which will be of great value for guiding future researches on VWA proteins and their important biological significance.

Ligand Binding Mechanism and Its Relationship with Conformational Changes in Adenine Riboswitch

Riboswitches are naturally occurring RNA aptamers that control the expression of essential bacterial genes by binding to specific small molecules. The binding with both high affinity and specificity induces conformational changes. Thus, riboswitches were proposed as a possible molecular target for developing antibiotics and chemical tools. The adenine riboswitch can bind not only to purine analogues but also to pyrimidine analogues. Here, long molecular dynamics (MD) simulations and molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) computational methodologies were carried out to show the differences in the binding model and the conformational changes upon five ligands binding. The binding free energies of the guanine riboswitch aptamer with C74U mutation complexes were compared to the binding free energies of the adenine riboswitch (AR) aptamer complexes. The calculated results are in agreement with the experimental data. The differences for the same ligand binding to two different aptamers are related to the electrostatic contribution. Binding dynamical analysis suggests a flexible binding pocket for the pyrimidine ligand in comparison with the purine ligand. The 18 μs of MD simulations in total indicate that both ligand-unbound and ligand-bound aptamers transfer their conformation between open and closed states. The ligand binding obviously affects the conformational change. The conformational states of the aptamer are associated with the distance between the mass center of two key nucleotides (U51 and A52) and the mass center of the other two key nucleotides (C74 and C75). The results suggest that the dynamical character of the binding pocket would affect its biofunction. To design new ligands of the adenine riboswitch, it is recommended to consider the binding affinities of the ligand and the conformational change of the ligand binding pocket.

Consistent Prediction of Mutation Effect on Drug Binding in HIV-1 Protease Using Alchemical Calculations

Despite a large number of antiretroviral drugs targeting HIV-1 protease for inhibition, mutations in this protein during the course of patient treatment can render them inefficient. This emerging resistance inspired numerous computational studies of the HIV-1 protease aimed at predicting the effect of mutations on drug binding in terms of free binding energy Δ G, as well as in mechanistic terms. In this study, we analyze ten different protease-inhibitor complexes carrying major resistance-associated mutations (RAMs) G48V, I50V, and L90M using molecular dynamics simulations. We demonstrate that alchemical free energy calculations can consistently predict the effect of mutations on drug binding. By explicitly probing different protonation states of the catalytic aspartic dyad, we reveal the importance of the correct choice of protonation state for the accuracy of the result. We also provide insight into how different mutations affect drug binding in their specific ways, with the unifying theme of how all of them affect the crucial drug binding regions of the protease.

Exploring the crizotinib resistance mechanism of NSCLC with the L1196M mutation using molecular dynamics simulation

Crizotinib is an anticancer tyrosine kinase inhibitor that is approved for use as a first-line treatment for some non-small-cell lung cancers. L1196M is the most frequently observed mutation in NSCLC patients. This mutation, known as the gatekeeper mutation in the ALK kinase domain, confers resistance to crizotinib by sterically blocking the binding of the drug. However, the molecular mechanism of crizotinib resistance caused by the L1196M mutation is still unclear. Molecular dynamics simulation was therefore utilized in this study to investigate the mechanism by which the L1196M mutation may affect crizotinib resistance. Our results suggest that larger fluctuations in some important regions of the mutant complex compared to the wild-type complex may contribute to the resistance of the mutant complex to crizotinib. Also, mutation-induced alterations to the secondary structure of the complex as well as unstable hydrogen-bonding patterns in the A-loop and P-loop regions decrease the total binding energy of the complex. This study therefore provides a molecular explanation for the resistance to crizotinib caused by the L1196M mutation, which could aid the design of more efficient and selective drugs.

Ligand Selectivity Mechanism and Conformational Changes in Guanine Riboswitch by Molecular Dynamics Simulations and Free Energy Calculations

Riboswitches regulate gene expression through direct and specific interactions with small metabolite molecules. Binding of a ligand to its RNA target is high selectivity and affinity and induces conformational changes of the RNA's secondary and tertiary structure. The structural difference of two purine riboswitches aptamers is caused by only one single mutation, where cytosine 74 in the guanine riboswitch is corresponding to a uracil 74 in adenine riboswitch. Here we employed molecular dynamics (MD) simulation, molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) and thermodynamic integration computational methodologies to evaluate the energetic and conformational changes of ligands binding to purine riboswitches. The snapshots used in MM-PBSA calculation were extracted from ten 50 ns MD simulation trajectories for each complex. These free energy results are in consistent with the experimental data and rationalize the selectivity of the riboswitches for different ligands. In particular, it is found that the loss in binding free energy upon mutation is mainly electrostatic in guanine (GUA) and riboswitch complex. Furthermore, new hydrogen bonds are found in mutated complexes. To reveal the conformational properties of guanine riboswitch, we performed a total of 6 μs MD simulations in both the presence and the absence of the ligand GUA. The MD simulations suggest that the conformation of guanine riboswitch depends on the distance of two groups in the binding pocket of ligand. The conformation is in a close conformation when U51-A52 is close to C74-U75.

Molecular dynamics and MM/GBSA-integrated protocol probing the correlation between biological activities and binding free energies of HIV-1 TAR RNA inhibitors

The interaction of HIV-1 transactivator protein Tat with its cognate transactivation response (TAR) RNA has emerged as a promising target for developing antiviral compounds and treating HIV infection, since it is a crucial step for efficient transcription and replication. In the present study, molecular dynamics (MD) simulations and MM/GBSA calculations have been performed on a series of neamine derivatives in order to estimate appropriate MD simulation time for acceptable correlation between ΔG_bind and experimental pIC₅₀ values. Initially, all inhibitors were docked into the active site of HIV-1 TAR RNA. Later to explore various conformations and examine the docking results, MD simulations were carried out. Finally, binding free energies were calculated using MM/GBSA method and were correlated with experimental pIC₅₀ values at different time scales (0-1 to 0-10 ns). From this study, it is clear that in case of neamine derivatives as simulation time increased the correlation between binding free energy and experimental pIC₅₀ values increased correspondingly. Therefore, the binding energies which can be interpreted at longer simulation times can be used to predict the bioactivity of new neamine derivatives. Moreover, in this work, we have identified some plausible critical nucleotide interactions with neamine derivatives that are responsible for potent inhibitory activity. Furthermore, we also provide some insights into a new class of oxadiazole-based back bone cyclic peptides designed by incorporating the structural features of neamine derivatives. On the whole, this approach can provide a valuable guidance for designing new potent inhibitors and modify the existing compounds targeting HIV-1 TAR RNA.

The conformational feasibility for the formation of reaching dimer in ASV and HIV integrase: a molecular dynamics study

Retroviral integrases are reported to form alternate dimer assemblies like the core-core dimer and reaching dimer. The core-core dimer is stabilized predominantly by an extensive interface between two catalytic core domains. The reaching dimer is stabilized by N-terminal domains that reach to form intermolecular interfaces with the other subunit's core and C-terminal domains (CTD), as well as CTD-CTD interactions. In this study, molecular dynamics (MD), Brownian dynamics (BD) simulations, and free energy analyses, were performed to elucidate determinants for the stability of the reaching dimer forms of full-length Avian Sarcoma Virus (ASV) and Human Immunodeficiency Virus (HIV) IN, and to examine the role of the C-tails (the last ~16-18 residues at the C-termini) in their structural dynamics. The dynamics of an HIV reaching dimer derived from small angle X-ray scattering and protein crosslinking data, was compared with the dynamics of a core-core dimer model derived from combining the crystal structures of two-domain fragments. The results showed that the core domains in the ASV reaching dimer express free dynamics, whereas those in the HIV reaching dimer are highly stable. BD simulations suggest a higher rate of association for the HIV core-core dimer than the reaching dimer. The predicted stability of these dimers was therefore ranked in the following order: ASV reaching dimer < HIV reaching dimer < composite core-core dimer. Analyses of MD trajectories have suggested residues that are critical for intermolecular contacts in each reaching dimer. Tests of these predictions and insights gained from these analyses could reveal a potential pathway for the association and dissociation of full-length IN multimers.

Computational Studies of a Mechanism for Binding and Drug Resistance in the Wild Type and Four Mutations of HIV-1 Protease with a GRL-0519 Inhibitor

Drug resistance of mutations in HIV-1 protease (PR) is the most severe challenge to the long-term efficacy of HIV-1 PR inhibitor in highly active antiretroviral therapy. To elucidate the molecular mechanism of drug resistance associated with mutations (D30N, I50V, I54M, and V82A) and inhibitor (GRL-0519) complexes, we have performed five molecular dynamics (MD) simulations and calculated the binding free energies using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method. The ranking of calculated binding free energies is in accordance with the experimental data. The free energy spectra of each residue and inhibitor interaction for all complexes show a similar binding model. Analysis based on the MD trajectories and contribution of each residues show that groups R2 and R3 mainly contribute van der Waals energies, while groups R1 and R4 contribute electrostatic interaction by hydrogen bonds. The drug resistance of D30N can be attributed to the decline in binding affinity of residues 28 and 29. The size of Val50 is smaller than Ile50 causes the residue to move, especially in chain A. The stable hydrophobic core, including the side chain of Ile54 in the wild type (WT) complex, became unstable in I54M because the side chain of Met54 is flexible with two alternative conformations. The binding affinity of Ala82 in V82A decreases relative to Val82 in WT. The present study could provide important guidance for the design of a potent new drug resisting the mutation inhibitors.

Predicting the impact of mutations on the specific activity of Bacillus thermocatenulatus lipase using a combined approach of docking and molecular dynamics

Lipases are important biocatalysts owing to their ability to catalyze diverse reactions with exceptional substrate specificities. A combined docking and molecular dynamics (MD) approach was applied to study the chain-length selectivity of Bacillus thermocatenulatus lipase (BTL2) towards its natural substrates (triacylglycerols). A scoring function including electrostatic, van der Waals (vdW) and desolvation energies along with conformational entropy was developed to predict the impact of mutation. The native BTL2 and its 6 mutants (F17A, V175A, V175F, D176F, T178V and I320F) were experimentally analyzed to determine their specific activities towards tributyrin (C4) or tricaprylin (C8), which were used to test our approach. Our scoring methodology predicted the chain-length selectivity of BTL2 with 85.7% (6/7) accuracy with a positive correlation between the calculated scores and the experimental activity values (r = 0.82, p = 0.0004). Additionally, the impact of mutation on activity was predicted with 75% (9/12) accuracy. The described study represents a fast and reliable approach to accurately predict the effect of mutations on the activity and selectivity of lipases and also of other enzymes. Copyright © 2016 John Wiley & Sons, Ltd.

E-pharmacophore filtering and molecular dynamics simulation studies in the discovery of potent drug-like molecules for chronic kidney disease

Chronic kidney disease (CKD) is a prominent health issue reported globally. The level of the vitamin D receptor (VDR) and cytochrome P450 enzyme 24-hydroxylase (CYP24A1) are crucial in the pathogenesis of secondary hyperparathyroidism (sHPT) in CKD. An elevated expression of the CYP24A1 leads to the deficiency of vitamin D and resistance to vitamin D therapy. Hence, VDR agonists and CYP24A1 antagonists are suggested to CKD patients for the management of biochemical complications. CTA-018 is a recently reported analog and acts as a potent CYP24A1 inhibitor. It inhibits CYP24A1 with an IC₅₀ 27 ± 6 nM, about 10 times more potentially than the non-selective inhibitor ketoconazole (253 ± 20 nM), and it is also been reported to induce the VDR expression. Thus, CTA-018 is under clinical trial among CKD patients. In this study, combined molecular docking and pharmacophore filtering were employed to identify compounds better than CTA-018. A huge set of 9127 compounds from Sweet Lead database were docked into the active site of VDR using Glide XP program. E-pharmacophore was developed from both the targets along with CTA-018. The compounds retrieved from the two different pharmacophore-based screening were re-docked into the active site of CYP24A1. The hits that bind well at both the active sites and matched with the pharmacophore models were considered as possible dual functional molecules against VDR and CYP24A1. Further, molecular dynamics simulation and subsequent energy decomposition analyses were also performed to study the role of specific amino acids in the active site of both VDR and CYP24A1.

Binding Free Energy Calculations of Nine FDA-approved Protease Inhibitors Against HIV-1 Subtype C I36T↑T Containing 100 Amino Acids Per Monomer

In this work, have investigated the binding affinities of nine FDA-approved protease inhibitor drugs against a new HIV-1 subtype C mutated protease, I36T↑T. Without an X-ray crystal structure, homology modelling was used to generate a three-dimensional model of the protease. This and the inhibitor models were employed to generate the inhibitor/I36T↑T complexes, with the relative positions of the inhibitors being superimposed and aligned using the X-ray crystal structures of the inhibitors/HIV-1 subtype B complexes as a reference. Molecular dynamics simulations were carried out on the complexes to calculate the average binding free energies for each inhibitor using the molecular mechanics generalized Born surface area (MM-GBSA) method. When compared to the binding free energies of the HIV-1 subtype B and subtype C proteases (calculated previously by our group using the same method), it was clear that the I36T↑T proteases mutations and insertion had a significant negative effect on the binding energies of the non-pepditic inhibitors nelfinavir, darunavir and tipranavir. On the other hand, ritonavir, amprenavir and indinavir show improved calculated binding energies in comparison with the corresponding data for wild-type C-SA protease. The computational model used in this study can be used to investigate new mutations of the HIV protease and help in establishing effective HIV drug regimes and may also aid in future protease drug design.

Water Sites, Networks, And Free Energies with Grand Canonical Monte Carlo

Water molecules play integral roles in the formation of many protein-ligand complexes, and recent computational efforts have been focused on predicting the thermodynamic properties of individual waters and how they may be exploited in rational drug design. However, when water molecules form highly coupled hydrogen-bonding networks, there is, as yet, no method that can rigorously calculate the free energy to bind the entire network or assess the degree of cooperativity between waters. In this work, we report theoretical and methodological developments to the grand canonical Monte Carlo simulation technique. Central to our results is a rigorous equation that can be used to calculate efficiently the binding free energies of water networks of arbitrary size and complexity. Using a single set of simulations, our methods can locate waters, estimate their binding affinities, capture the cooperativity of the water network, and evaluate the hydration free energy of entire protein binding sites. Our techniques have been applied to multiple test systems and compare favorably to thermodynamic integration simulations and experimental data. The implications of these methods in drug design are discussed.

Revealing the binding modes and the unbinding of 14-3-3σ proteins and inhibitors by computational methods

The 14-3-3σ proteins are a family of ubiquitous conserved eukaryotic regulatory molecules involved in the regulation of mitogenic signal transduction, apoptotic cell death, and cell cycle control. A lot of small-molecule inhibitors have been identified for 14-3-3 protein-protein interactions (PPIs). In this work, we carried out molecular dynamics (MD) simulations combined with molecular mechanics generalized Born surface area (MM-GBSA) method to study the binding mechanism between a 14-3-3σ protein and its eight inhibitors. The ranking order of our calculated binding free energies is in agreement with the experimental results. We found that the binding free energies are mainly from interactions between the phosphate group of the inhibitors and the hydrophilic residues. To improve the binding free energy of Rx group, we designed the inhibitor R9 with group R9 = 4-hydroxypheny. However, we also found that the binding free energy of inhibitor R9 is smaller than that of inhibitor R1. By further using the steer molecular dynamics (SMD) simulations, we identified a new hydrogen bond between the inhibitor R8 and residue Arg64 in the pulling paths. The information obtained from this study may be valuable for future rational design of novel inhibitors, and provide better structural understanding of inhibitor binding to 14-3-3σ proteins.

Searching the Sequence Space for Potent Aptamers Using SELEX in Silico

To isolate functional nucleic acids that bind to defined targets with high affinity and specificity, which are known as aptamers, the systematic evolution of ligands by exponential enrichment (SELEX) methodology has emerged as the preferred approach. Here, we propose a computational approach, SELEX in silico, that allows the sequence space to be more thoroughly explored regarding binding of a certain target. Our approach consists of two steps: (i) secondary structure-based sequence screening, which aims to collect the sequences that can form a desired RNA motif as an enhanced initial library, followed by (ii) sequence enrichment regarding target binding by molecular dynamics simulation-based virtual screening. Our SELEX in silico method provided a practical computational solution to three key problems in aptamer sequence searching: design of nucleic acid libraries, knowledge of sequence enrichment, and identification of potent aptamers. Six potent theophylline-binding aptamers, which were isolated by SELEX in silico from a sequence space containing 4(13) sequences, were experimentally verified to bind theophylline with high affinity: Kd ranging from 0.16 to 0.52 μM, compared with the dissociation constant of the original aptamer-theophylline, 0.32 μM. These results demonstrate the significant potential of SELEX in silico as a new method for aptamer discovery and optimization.

Molecular dynamics study of biodegradation of azo dyes via their interactions with AzrC azoreductase

Azo dyes are one of the most important class of dyes, which have been widely used in industries. Because of the environmental pollution of azo dyes, many studies have been performed to study their biodegradation using bacterial systems. In present work, the AzrC of mesophilic gram-positive Bacillus sp. B29 has been considered to study its interaction with five common azo dyes (orange G, acid red 88, Sudan I, orange I, and methyl red). The molecular dynamics simulations have been employed to study the interaction between AzrC and azo dyes. The trajectory was confirmed using root mean square deviation and the root mean square fluctuation analyses. Then, the hydrogen bond and alanine scanning analyses were performed to reveal active site residues. Phe105 (A), Phe125 (B), Phe172 (B), and Pro132 (B) have been found as the most important hydrophobic residues whereas Asn104 (A), Tyr127 (B), and Asn187 (A) have key role in making hydrogen bond. The results of molecular mechanics Poisson-Boltzmann surface area and molecular mechanics generalized Born surface area calculations proved that the hydrophobic azo dyes like Acid red 88 binds more tightly to the AzrC protein. The calculated data suggested MR A 121 (B) I as a potential candidate for improving the AzrC-MR interactions.

Biological and computational evaluation of resveratrol inhibitors against Alzheimer's disease

It has been reported that beta amyloid induces production of radical oxygen species and oxidative stress in neuronal cells, which in turn upregulates β-secretase (BACE-1) expression and beta amyloid levels, thereby propagating oxidative stress and increasing neuronal injury. A series of resveratrol derivatives, known to be inhibitors of oxidative stress-induced neuronal cell death (oxytosis) were biologically evaluated against BACE-1 using homogeneous time-resolved fluorescence (TRF) assay. Correlation between oxytosis inhibitory and BACE-1 inhibitory activity of resveratrol derivatives was statistically significant, supporting the notion that BACE-1 may act as pivotal mediator of neuronal cell oxytosis. Four of the biologically evaluated resveratrol analogs demonstrated considerably higher activity than resveratrol in either assay. The discovery of some "hits" led us to initiate detailed docking studies associated with Molecular Dynamics in order to provide a plausible explanation for the experimental results and understand their molecular basis of action.

Rapid, Precise, and Reproducible Prediction of Peptide-MHC Binding Affinities from Molecular Dynamics That Correlate Well with Experiment

The presentation of potentially pathogenic peptides by major histocompatibility complex (MHC) molecules is one of the most important processes in adaptive immune defense. Prediction of peptide-MHC (pMHC) binding affinities is therefore a principal objective of theoretical immunology. Machine learning techniques achieve good results if substantial experimental training data are available. Approaches based on structural information become necessary if sufficiently similar training data are unavailable for a specific MHC allele, although they have often been deemed to lack accuracy. In this study, we use a free energy method to rank the binding affinities of 12 diverse peptides bound by a class I MHC molecule HLA-A*02:01. The method is based on enhanced sampling of molecular dynamics calculations in combination with a continuum solvent approximation and includes estimates of the configurational entropy based on either a one or a three trajectory protocol. It produces precise and reproducible free energy estimates which correlate well with experimental measurements. If the results are combined with an amino acid hydrophobicity scale, then an extremely good ranking of peptide binding affinities emerges. Our approach is rapid, robust, and applicable to a wide range of ligand-receptor interactions without further adjustment.

Pentacycloundecane lactam vs lactone norstatine type protease HIV inhibitors: binding energy calculations and DFT study

Background

Novel pentacycloundecane (PCU)-lactone-CO-EAIS peptide inhibitors were designed, synthesized, and evaluated against wild-type C-South African (C-SA) HIV-1 protease. Three compounds are reported herein, two of which displayed IC₅₀ values of less than 1.00 μM. A comparative MM-PB(GB)SA binding free energy of solvation values of PCU-lactam and lactone models and their enantiomers as well as the PCU-lactam-NH-EAIS and lactone-CO-EAIS peptide inhibitors and their corresponding diastereomers complexed with South African HIV protease (C-SA) was performed. This will enable us to rationalize the considerable difference between inhibitory concentration (IC₅₀) of PCU-lactam-NH-EAIS and PCU-lactone-CO-EAIS peptides.

Results

The PCU-lactam model exhibited more negative calculated binding free energies of solvation than the PCU-lactone model. The same trend was observed for the PCU-peptide inhibitors, which correspond to the experimental activities for the PCU-lactam-NH-EAIS peptide (IC₅₀ = 0.076 μM) and the PCU-lactone-CO-EAIS peptide inhibitors (IC₅₀ = 0.850 μM). Furthermore, a density functional theory (DFT) study on the natural atomic charges of the nitrogen and oxygen atoms of the three PCU-lactam, PCU-lactim and PCU-lactone models were performed using natural bond orbital (NBO) analysis. Electrostatic potential maps were also used to visualize the electron density around electron-rich regions. The asymmetry parameter (η) and quadrupole coupling constant (χ) values of the nitrogen and oxygen nuclei of the model compounds were calculated at the same level of theory. Electronic molecular properties including polarizability and electric dipole moments were also calculated and compared. The Gibbs theoretical free solvation energies of solvation (∆G_solv) were also considered.

Conclusions

A general trend is observed that the lactam species appears to have a larger negative charge distribution around the heteroatoms, larger quadrupole constant, dipole moment and better solvation energy, in comparison to the PCU-lactone model. It can be argued that these characteristics will ensure better eletronic interaction between the lactam and the receptor, corresponding to the observed HIV protease activities in terms of experimental IC₅₀ data.

Electronic supplementary material

The online version of this article (doi:10.1186/s12929-015-0115-5) contains supplementary material, which is available to authorized users.

The impact of active site mutations of South African HIV PR on drug resistance: Insight from molecular dynamics simulations, binding free energy and per-residue footprints

Molecular dynamics simulations and binding free energy calculations were used to provide an understanding of the impact of active site drug-resistant mutations of the South African HIV protease subtype C (C-SA HIV PR), V82A and V82F/I84V on drug resistance. Unique per-residue interaction energy 'footprints' were developed to map the overall drug-binding profiles for the wild type and mutants. Results confirmed that these mutations altered the overall binding landscape of the amino acid residues not only in the active site region but also in the flaps as well. Four FDA-approved drugs were investigated in this study; these include ritonavir (RTV), saquinavir (SQV), indinavir (IDV), and nelfinavir (NFV). Computational results compared against experimental findings were found to be complementary. Against the V82F/I84V variant, saquinavir, indinavir, and nelfinavir lose remarkable entropic contributions relative to both wild-type and V82A C-SA HIV PRs. The per-residue energy 'footprints' and the analysis of ligand-receptor interactions for the drug complexes with the wild type and mutants have also highlighted the nature of drug interactions. The data presented in this study will prove useful in the design of more potent inhibitors effective against drug-resistant HIV strains.

Target specific proteochemometric model development for BACE1 – protein flexibility and structural water are critical in virtual screening

Structural water and protein plasticity are important factors for BACE1 targeted ligand virtual screening.

Exploring the drug resistance of V32I and M46L mutant HIV-1 protease to inhibitor TMC114: flap dynamics and binding mechanism

Inhibitors of HIV-1 protease (HIV-1-pr) generally only bind to the active site of the protease. However, for some mutants such as V32I and M46L the TMC114 can bind not only to the active cavity but also to the groove of the flexible flaps. Although the second binding site suggests the higher efficiency of the drug against HIV-1-pr, the drug resistance in HIV-1-pr due to mutations cannot be ignored, which prompts us to investigate the molecular mechanisms of drug resistance and behavior of double bound TMC114 (2T) to HIV-1-pr. The conformational dynamics of HIV-1-pr and the binding of TMC114 to the WT, V32I and M46L mutants were investigated with all-atom molecular dynamic (MD) simulation. The 20 ns MD simulation shows many fascinating effects of the inhibitor binding to the WT and mutant proteases. MM-PBSA calculations explain the binding free energies unfavorable for the M46L and V32I mutants as compared to the WT. For the single binding (1T) the less binding affinity can be attributed to the entropic loss for both V32I-1T and M46L-1T. Although the second binding of TMC114 with flap does increase binding energy for the mutants (V32I-2T and M46L-2T), the considerable entropy loss results in the lower binding Gibbs free energies. Thus, binding of TMC114 in the flap region does not help much in the total gain in binding affinity of the system, which was verified from this study and thereby validating experiments.

Impact of M36I polymorphism on the interaction of HIV-1 protease with its substrates: insights from molecular dynamics

Background

Over the last decades, a vast structural knowledge has been gathered on the HIV-1 protease (PR). Noticeably, most of the studies focused the B-subtype, which has the highest prevalence in developed countries. Accordingly, currently available anti-HIV drugs target this subtype, with considerable benefits for the corresponding patients.

However, in developing countries, there is a wide variety of HIV-1 subtypes carrying PR polymorphisms related to reduced drug susceptibility. The non-active site mutation, M36I, is the most frequent polymorphism, and is considered as a non-B subtype marker.

Yet, the structural impact of this substitution on the PR structure and on the interaction with natural substrates remains poorly documented.

Results

Herein, we used molecular dynamics simulations to investigate the role of this polymorphism on the interaction of PR with six of its natural cleavage-sites substrates.

Free energy analyses by MMPB/SA calculations showed an affinity decrease of M36I-PR for the majority of its substrates. The only exceptions were the RT-RH, with equivalent affinity, and the RH-IN, for which an increased affinity was found. Furthermore, molecular simulations suggest that, unlike other peptides, RH-IN induced larger structural fluctuations in the wild-type enzyme than in the M36I variant.

Conclusions

With multiple approaches and analyses we identified structural and dynamical determinants associated with the changes found in the binding affinity of the M36I variant. This mutation influences the flexibility of both PR and its complexed substrate. The observed impact of M36I, suggest that combination with other non-B subtype polymorphisms, could lead to major effects on the interaction with the 12 known cleavage sites, which should impact the virion maturation.

Systematic molecular dynamics, MM-PBSA, and ab initio approaches to the saquinavir resistance mechanism in HIV-1 PR due to 11 double and multiple mutations

Mutations in the human immunodeficiency virus (HIV) enable virus replication even when appropriate antiretroviral therapy is followed, thus leading to the emergence of drug resistance. In a previous work, we systematically examined seven single mutations that are associated with saquinavir (SQV) resistance in HIV-1 protease (Tzoupis, H.; Leonis, G.; Mavromoustakos, T.; Papadopoulos, M. G. J. Chem. Theory Comput. 2013, 9, 1754-1764). Herein, we extend our analysis, which includes seven double (G48V-V82A, L10I-G48V, G48V-L90M, I84V-L90M, L10I-V82A, L10I-L63P, A71V-G73S) and four multiple (L10I-L63P-A71V, L10I-G48V-V82A, G73S-I84V-L90M, L10I-L63P-A71V-G73S-I84V-L90M) SQV-HIV-1 PR mutant complexes, in an attempt to generalize our findings and formulate the main elements of the SQV resistance mechanism in the protease. On the basis of molecular dynamics (MD), molecular mechanics Poisson-Boltzmann surface area (MM-PBSA), and ab initio computational approaches, we identified specific features that constitute the HIV-1 PR mechanism of resistance at the molecular level: the low flexibility of SQV in the binding cavity and the preservation of hydrogen bonding (HB) and van der Waals interactions between SQV and several active-site (Gly27/27', Asp29/29'/30/30', especially Asp25/25') and flap (Ile50/50', Gly48/48') residues of the protease contribute significantly to efficient binding. The total enthalpy loss in all mutants is mostly due to the loss in enthalpy of the active-site region. Furthermore, it was observed that mutation accumulation may induce stabilization to SQV and to the flaps through enhanced HB interactions that lead to improved inhibition (e.g., accumulation of mutations in complexes containing L10I, G48V, L63P, I84V, or L90M single mutations). It was also concluded that permanent flap closure is obtained independently of mutations and SQV binding is mostly driven by van der Waals, nonpolar, and exchange-energy contributions. Importantly, it was indicated that the optimal positioning of SQV and the structure of the binding cavity are tightly coupled, since small changes in geometry may affect the binding energy greatly. The results of our theoretical approaches are in agreement with experimental evidence and provide a reliable description of SQV resistance in HIV-1 PR.

Interaction investigations of HipA binding to HipB dimer and HipB dimer + DNA complex: a molecular dynamics simulation study

We carried out molecular dynamics simulations and free energy calculations for a series of ternary and diplex models for the HipA protein, HipB dimer, and DNA molecule to address the mechanism of HipA sequestration and the binding order of events from apo HipB/HipA to 2HipA + HipB dimer + DNA complex. The results revealed that the combination of DNA with the HipB dimer is energetically favorable for the combination of HipB dimer with HipA protein. The binding of DNA to HipB dimer induces a long-range allosteric communication from the HipB2 -DNA interface to the HipA-HipB2 interface, which involves the closeness of α1 helices of HipB dimer to HipA protein and formations of extra hydrogen bonds in the HipA-HipB2 interface through the extension of α2/3 helices in the HipB dimer. These simulated results suggested that the DNA molecule, as a regulative media, modulates the HipB dimer conformation, consequently increasing the interactions of HipB dimer with the HipA proteins, which explains the mechanism of HipA sequestration reported by the previous experiment. Simultaneously, these simulations also explored that the thermodynamic binding order in a simulated physiological environment, that is, the HipB dimer first bind to DNA to form HipB dimer + DNA complex, then capturing strongly the HipA proteins to form a ternary complex, 2HipA + HipB dimer + DNA, for sequestrating HipA in the nucleoid.

Computing Clinically Relevant Binding Free Energies of HIV-1 Protease Inhibitors

The use of molecular simulation to estimate the strength of macromolecular binding free energies is becoming increasingly widespread, with goals ranging from lead optimization and enrichment in drug discovery to personalizing or stratifying treatment regimes. In order to realize the potential of such approaches to predict new results, not merely to explain previous experimental findings, it is necessary that the methods used are reliable and accurate, and that their limitations are thoroughly understood. However, the computational cost of atomistic simulation techniques such as molecular dynamics (MD) has meant that until recently little work has focused on validating and verifying the available free energy methodologies, with the consequence that many of the results published in the literature are not reproducible. Here, we present a detailed analysis of two of the most popular approximate methods for calculating binding free energies from molecular simulations, molecular mechanics Poisson-Boltzmann surface area (MMPBSA) and molecular mechanics generalized Born surface area (MMGBSA), applied to the nine FDA-approved HIV-1 protease inhibitors. Our results show that the values obtained from replica simulations of the same protease-drug complex, differing only in initially assigned atom velocities, can vary by as much as 10 kcal mol-1, which is greater than the difference between the best and worst binding inhibitors under investigation. Despite this, analysis of ensembles of simulations producing 50 trajectories of 4 ns duration leads to well converged free energy estimates. For seven inhibitors, we find that with correctly converged normal mode estimates of the configurational entropy, we can correctly distinguish inhibitors in agreement with experimental data for both the MMPBSA and MMGBSA methods and thus have the ability to rank the efficacy of binding of this selection of drugs to the protease (no account is made for free energy penalties associated with protein distortion leading to the over estimation of the binding strength of the two largest inhibitors ritonavir and atazanavir). We obtain improved rankings and estimates of the relative binding strengths of the drugs by using a novel combination of MMPBSA/MMGBSA with normal mode entropy estimates and the free energy of association calculated directly from simulation trajectories. Our work provides a thorough assessment of what is required to produce converged and hence reliable free energies for protein-ligand binding.

Molecular modeling and residue interaction network studies on the mechanism of binding and resistance of the HCV NS5B polymerase mutants to VX-222 and ANA598

Hepatitis C virus (HCV) NS5B protein is an RNA-dependent RNA polymerase (RdRp) with essential functions in viral genome replication and represents a promising therapeutic target to develop direct-acting antivirals (DAAs). Multiple nonnucleoside inhibitors (NNIs) binding sites have been identified within the polymerase. VX-222 and ANA598 are two NNIs targeting thumb II site and palm I site of HCV NS5B polymerase, respectively. These two molecules have been shown to be very effective in phase II clinical trials. However, the emergence of resistant HCV replicon variants (L419M, M423T, I482L mutants to VX-222 and M414T, M414L, G554D mutants to ANA598) has significantly decreased their efficacy. To elucidate the molecular mechanism about how these mutations influenced the drug binding mode and decreased drug efficacy, we studied the binding modes of VX-222 and ANA598 to wild-type and mutant polymerase by molecular modeling approach. Molecular dynamics (MD) simulations results combined with binding free energy calculations indicated that the mutations significantly altered the binding free energy and the interaction for the drugs to polymerase. The further per-residue binding free energy decomposition analysis revealed that the mutations decreased the interactions with several key residues, such as L419, M423, L474, S476, I482, L497, for VX-222 and L384, N411, M414, Y415, Q446, S556, G557 for ANA598. These were the major origins for the resistance to these two drugs. In addition, by analyzing the residue interaction network (RIN) of the complexes between the drugs with wild-type and the mutant polymerase, we found that the mutation residues in the networks involved in the drug resistance possessed a relatively lower size of topology centralities. The shift of betweenness and closeness values of binding site residues in the mutant polymerase is relevant to the mechanism of drug resistance of VX-222 and ANA598. These results can provide an atomic-level understanding about the mechanisms of drug resistance conferred by the studied mutations and will be helpful to design more potent inhibitors which could effectively overcome drug resistance of antivirus agents.

Ligand binding and dynamics of the monomeric epidermal growth factor receptor ectodomain

The ectodomain of the human epidermal growth factor receptor (hEGFR) controls input to several cell signalling networks via binding with extracellular growth factors. To gain insight into the dynamics and ligand binding of the ectodomain, the hEGFR monomer was subjected to molecular dynamics simulation. The monomer was found to be substantially more flexible than the ectodomain dimer studied previously. Simulations where the endogeneous ligand EGF binds to either Subdomain I or Subdomain III, or where hEGFR is unbound, show significant differences in dynamics. The molecular mechanics Poisson-Boltzmann surface area method has been used to derive relative free energies of ligand binding, and we find that the ligand is capable of binding either subdomain with a slight preference for III. Alanine-scanning calculations for the effect of selected ligand mutants on binding reproduce the trends of affinity measurements. Taken together, these results emphasize the possible role of the ectodomain monomer in the initial step of ligand binding, and add details to the static picture obtained from crystal structures.

Efficient Computation of Small-Molecule Configurational Binding Entropy and Free Energy Changes by Ensemble Enumeration

Here we present a novel, end-point method using the dead-end-elimination and A* algorithms to efficiently and accurately calculate the change in free energy, enthalpy, and configurational entropy of binding for ligand–receptor association reactions. We apply the new approach to the binding of a series of human immunodeficiency virus (HIV-1) protease inhibitors to examine the effect ensemble reranking has on relative accuracy as well as to evaluate the role of the absolute and relative ligand configurational entropy losses upon binding in affinity differences for structurally related inhibitors. Our results suggest that most thermodynamic parameters can be estimated using only a small fraction of the full configurational space, and we see significant improvement in relative accuracy when using an ensemble versus single-conformer approach to ligand ranking. We also find that using approximate metrics based on the single-conformation enthalpy differences between the global minimum energy configuration in the bound as well as unbound states also correlates well with experiment. Using a novel, additive entropy expansion based on conditional mutual information, we also analyze the source of ligand configurational entropy loss upon binding in terms of both uncoupled per degree of freedom losses as well as changes in coupling between inhibitor degrees of freedom. We estimate entropic free energy losses of approximately +24 kcal/mol, 12 kcal/mol of which stems from loss of translational and rotational entropy. Coupling effects contribute only a small fraction to the overall entropy change (1–2 kcal/mol) but suggest differences in how inhibitor dihedral angles couple to each other in the bound versus unbound states. The importance of accounting for flexibility in drug optimization and design is also discussed.

Computational study on the drug resistance mechanism against HCV NS3/4A protease inhibitors vaniprevir and MK-5172 by the combination use of molecular dynamics simulation, residue interaction network, and substrate envelope analysis

Hepatitis C virus (HCV) NS3/4A protease is an important and attractive target for anti-HCV drug development and discovery. Vaniprevir (phase III clinical trials) and MK-5172 (phase II clinical trials) are two potent antiviral compounds that target NS3/4A protease. However, the emergence of resistance to these two inhibitors reduced the effectiveness of vaniprevir and MK-5172 against viral replication. Among the drug resistance mutations, three single-site mutations at residues Arg155, Ala156, and Asp168 in NS3/4A protease are especially important due to their resistance to nearly all inhibitors in clinical development. A detailed understanding of drug resistance mechanism to vaniprevir and MK-5172 is therefore very crucial for the design of novel potent agents targeting viral variants. In this work, molecular dynamics (MD) simulation, binding free energy calculation, free energy decomposition, residue interaction network (RIN), and substrate envelope analysis were used to study the detailed drug resistance mechanism of the three mutants R155K, A156T, and D168A to vaniprevir and MK-5172. MD simulation was used to investigate the binding mode for these two inhibitors to wild-type and resistant mutants of HCV NS3/4A protease. Binding free energy calculation and free energy decomposition analysis reveal that drug resistance mutations reduced the interactions between the active site residues and substituent in the P2 to P4 linker of vaniprevir and MK-5172. Furthermore, RIN and substrate envelope analysis indicate that the studied mutations of the residues are located outside the substrate (4B5A) binding site and selectively decrease the affinity of inhibitors but not the activity of the enzyme and consequently help NS3/4A protease escape from the effect of the inhibitors without influencing the affinity of substrate binding. These findings can provide useful information for understanding the drug resistance mechanism against vaniprevir and MK-5172. The results can also provide some potential clues for further design of novel inhibitors that are less susceptible to drug resistance.

Patient-specific modelling in drug design, development and selection including its role in clinical decision-making

Genomics has made enormous progress in the twelve years since the publication of the first draft human genome sequence, but it has not yet been translated into the clinic. Despite spiralling development costs, the number of new drug registrations is not increasing. One reason for this lies in the genetic complexity of disease. Most diseases involve dysregulation in pathways that involve many genes, and many (including most cancers) are themselves genetically heterogeneous. Systems biology involves the multi-level simulation of physiology, cell biology and biochemistry using complex computational techniques. We show here using case studies in cancer and HIV how such computational models, and particularly models based on individual patient data, can be used for drug design and development, and in the selection of the appropriate treatment for a given patient in the face of resistance mutations. If these techniques are to be adopted in routine clinical practice, clinicians will need better training in modern approaches to the integrated analysis of large-scale heterogeneous data and multi-scale models, while developers will need to provide much more usable tools. Investment in computational infrastructure is needed so that results can be returned on clinically relevant timescales and data warehouses designed with data protection as well as accessibility in mind.

Importance of polar solvation and configurational entropy for design of antiretroviral drugs targeting HIV-1 protease

Both KNI-10033 and KNI-10075 are high affinity preclinical HIV-1 protease (PR) inhibitors with affinities in the picomolar range. In this work, the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method has been used to investigate the potency of these two HIV-1 PR inhibitors against the wild-type and mutated proteases assuming that potency correlates with the affinity of the drugs for the target protein. The decomposition of the binding free energy reveals the origin of binding affinities or mutation-induced affinity changes. Our calculations indicate that the mutation I50V causes drug resistance against both inhibitors. On the other hand, we predict that the mutant I84V causes drug resistance against KNI-10075 while KNI-10033 is more potent against the I84V mutant compared to wild-type protease. Drug resistance arises mainly from unfavorable shifts in van der Waals interactions and configurational entropy. The latter indicates that neglecting changes in configurational entropy in the computation of relative binding affinities as often done is not appropriate in general. For the bound complex PR(I50V)-KNI-10075, an increased polar solvation free energy also contributes to the drug resistance. The importance of polar solvation free energies is revealed when interactions governing the binding of KNI-10033 or KNI-10075 to the wild-type protease are compared to the inhibitors darunavir or GRL-06579A. Although the contributions from intermolecular electrostatic and van der Waals interactions as well as the nonpolar component of the solvation free energy are more favorable for PR-KNI-10033 or PR-KNI-10075 compared to PR-DRV or PR-GRL-06579A, both KNI-10033 and KNI-10075 show a similar affinity as darunavir and a lower binding affinity relative to GRL-06579A. This is because of the polar solvation free energy which is less unfavorable for darunavir or GRL-06579A relative to KNI-10033 or KNI-10075. The importance of the polar solvation as revealed here highlights that structural inspection alone is not sufficient for identifying the key contributions to binding affinities and affinity changes for the design of drugs but that solvation effects must be taken into account. A detailed understanding of the molecular forces governing binding and drug resistance might assist in the design of new inhibitors against HIV-1 PR variants that are resistant against current drugs.