Minimum sequence requirements for the binding of paromomycin to the rRNA decoding site A

We have recently introduced a computational methodology that combines molecular dynamics (MD) simulations, free-energy calculations, and in vitro binding assays to predict the minimum RNA structural requirements for selective, high-affinity RNA binding to small-molecule ligands. Here, we show that this methodology can be applied to the conformationally flexible aminoglycoside antibiotic paromomycin. A RNA consisting of an 11-mer:10-mer duplex that contains one 16S ribosome RNA decoding A-site bound to paromomycin was simulated for 4 ns. The methodology predicts that the 11-mer:10-mer duplex binds to paromomycin with high affinity, whereas smaller RNA duplexes lose complex stability and the ability to bind paromomycin. The predicted high-affinity binding to paromomycin of the 11-mer:10-mer duplex was confirmed experimentally (EC(50) = 0.28 microM), as well as the inability of smaller complexes to bind. Our simulations show good agreement with experiment for dynamic and structural properties of the isolated A-site, including hydrogen-bonding networks and RNA structural rearrangements upon ligand binding. The results suggest that MD simulations can supplement in vitro methods as a tool for predicting minimum RNA-binding motifs for both small, rigid ligands, and large, flexible ligands when structural information is available.     

Insights into the functional role of protonation states in the HIV-1 protease-BEA369 complex: molecular dynamics simulations and free energy calculations

The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method combined with molecular dynamics (MD) simulations were used to investigate the functional role of protonation in human immunodeficiency virus type 1 (HIV-1) protease complexed with the inhibitor BEA369. Our results demonstrate that protonation of two aspartic acids (Asp25/Asp25′) has a strong influence on the dynamics behavior of the complex, the binding free energy of BEA369, and inhibitor–residue interactions. Relative binding free energies calculated using the MM-PBSA method show that protonation of Asp25 results in the strongest binding of BEA369 to HIV-1 protease. Inhibitor–residue interactions computed by the theory of free energy decomposition also indicate that protonation of Asp25 has the most favorable effect on binding of BEA369. In addition, hydrogen-bond analysis based on the trajectories of the MD simulations shows that protonation of Asp25 strongly influences the water-mediated link of a conserved water molecule, Wat301. We expect that the results of this study will contribute significantly to binding calculations for BEA369, and to the design of high affinity inhibitors.     

Molecular dynamics simulations and MM–PBSA calculations of the lectin from snowdrop (Galanthus nivalis)

Galanthus nivalis agglutinin (GNA), a mannose-specific lectin from snowdrop bulbs, is a member of the monocot mannose-specific lectin family and exhibits antiviral activity toward HIV. In the present study, molecular dynamics (MD) simulations were performed to study the interaction between GNA and its carbohydrate ligand over a specific time span. By analysis of the secondary structures, it was observed that the GNA conformation maintains rather stable along the trajectories and the high fluctuations were only centered on the carbohydrate recognition domains. Our MD simulations also reproduced most of the hydrogen bonds observed in the x-ray crystal structure. Furthermore, the obtained MD trajectories were used to estimate the binding free energy of the complex using the molecular mechanics/Poisson Boltzmann surface area (MM-PBSA) method. It was revealed by the inspection of the binding free energy components that the major contributions to the complex stability arose from electrostatic interactions.     

Hemolytic mechanism of dioscin proposed by molecular dynamics simulations

Saponins are a class of compounds containing a triterpenoid or steroid core with some attached carbohydrate modules. Many saponins cause hemolysis. However, the hemolytic mechanism of saponins at the molecular level is not yet fully understood. In an attempt to explore this issue, we have studied dioscin—a saponin with high hemolytic activity—through extensive molecular dynamics (MD) simulations. Firstly, all-atom MD simulations of 8 ns duration were conducted to study the stability of the dioscin–cholesterol complex and the cholesterol–cholesterol complex in water and in decane, respectively. MM-GB/SA computations indicate that the dioscin–cholesterol complex is energetically more favorable than the cholesterol–cholesterol complex in a non-polar environment. Next, several coarse-grained MD simulations of 400 ns duration were conducted to directly observe the distribution of multiple dioscin molecules on a DPPC-POPC-PSM-CHOL lipid bilayer. Our results indicate that dioscin can penetrate into the lipid bilayer, accumulate in the lipid raft micro-domain, and then bind cholesterol. This leads to the destabilization of lipid raft and consequent membrane curvature, which may eventually result in the hemolysis of red cells. This possible mechanism of hemolysis can well explain some experimental observations on hemolysis.     

Markov State Models Reveal a Two-Step Mechanism of miRNA Loading into the Human Argonaute Protein: Selective Binding followed by Structural Re-arrangement

Argonaute (Ago) proteins and microRNAs (miRNAs) are central components in RNA interference, which is a key cellular mechanism for sequence-specific gene silencing. Despite intensive studies, molecular mechanisms of how Ago recognizes miRNA remain largely elusive. In this study, we propose a two-step mechanism for this molecular recognition: selective binding followed by structural re-arrangement. Our model is based on the results of a combination of Markov State Models (MSMs), large-scale protein-RNA docking, and molecular dynamics (MD) simulations. Using MSMs, we identify an open state of apo human Ago-2 in fast equilibrium with partially open and closed states. Conformations in this open state are distinguished by their largely exposed binding grooves that can geometrically accommodate miRNA as indicated in our protein-RNA docking studies. miRNA may then selectively bind to these open conformations. Upon the initial binding, the complex may perform further structural re-arrangement as shown in our MD simulations and eventually reach the stable binary complex structure. Our results provide novel insights in Ago-miRNA recognition mechanisms and our methodology holds great potential to be widely applied in the studies of other important molecular recognition systems.     

Advanced replica-exchange sampling to study the flexibility and plasticity of peptides and proteins

Molecular dynamics (MD) simulations are ideally suited to investigate protein and peptide plasticity and flexibility simultaneously at high spatial (atomic) and high time resolution. However, the applicability is still limited by the force field accuracy and by the maximum simulation time that can be routinely achieved in current MD simulations. In order to improve the sampling the replica-exchange (REMD) methodology has become popular and is now the most widely applied advanced sampling approach. Many variants of the REMD method have been designed to reduce the computational demand or to enhance sampling along specific sets of conformational variables. An overview on recent methodological advances and discussion of specific aims and advantages of the approaches will be given. Applications in the area of free energy simulations and advanced sampling of intrinsically disordered peptides and proteins will also be discussed. This article is part of a Special Issue entitled: The emerging dynamic view of proteins: Protein plasticity in allostery, evolution and self-assembly.     

Docking study and free energy simulation of the complex between p53 DNA-binding domain and azurin

Molecular interaction between p53 tumor suppressor and the copper protein azurin (AZ) has been demonstrated to enhance p53 stability and hence antitumoral function, opening new perspectives in cancer treatment. While some experimental work has provided evidence for AZ binding to p53, no crystal structure for the p53-AZ complex was solved thus far. In this work the association between AZ and the p53 DNA-binding domain (DBD) was investigated by computational methods. Using a combination of rigid-body protein docking, experimental mutagenesis information, and cluster analysis 10 main p53 DBD-AZ binding modes were generated. The resulting structures were further characterized by molecular dynamics (MD) simulations and free energy calculations. We found that the highest scored docking conformation for the p53 DBD-AZ complex also yielded the most favorable free energy value. This best three-dimensional model for the complex was validated by using a computational mutagenesis strategy. In this structure AZ binds to the flexible L(1) and s(7)-s(8) loops of the p53 DBD and stabilizes them through protein-protein tight packing interactions, resulting in high degree of both surface matching and electrostatic complementarity.     

Multi‐scale characterization of the energy landscape of proteins with application to the C3D/Efb‐C complex

We present a novel multi‐level methodology to explore and characterize the low energy landscape and the thermodynamics of proteins. Traditional conformational search methods typically explore only a small portion of the conformational space of proteins and are hard to apply to large proteins due to the large amount of calculations required. In our multi‐scale approach, we first provide an initial characterization of the equilibrium state ensemble of a protein using an efficient computational conformational sampling method. We then enrich the obtained ensemble by performing short Molecular Dynamics (MD) simulations on selected conformations from the ensembles as starting points. To facilitate the analysis of the results, we project the resulting conformations on a low‐dimensional landscape to efficiently focus on important interactions and examine low energy regions. This methodology provides a more extensive sampling of the low energy landscape than an MD simulation starting from a single crystal structure as it explores multiple trajectories of the protein. This enables us to obtain a broader view of the dynamics of proteins and it can help in understanding complex binding, improving docking results and more. In this work, we apply the methodology to provide an extensive characterization of the bound complexes of the C3d fragment of human Complement component C3 and one of its powerful bacterial inhibitors, the inhibitory domain of Staphylococcus aureus extra‐cellular fibrinogen‐binding domain (Efb‐C) and two of its mutants. We characterize several important interactions along the binding interface and define low free energy regions in the three complexes. Proteins 2010. © 2009 Wiley‐Liss, Inc.     

Novel insights into the dynamics behavior of glucagon-like peptide-1 receptor with its small molecule agonists

Abstract

The glucagon-like peptide-1 receptor (GLP-1R) is a well-known target of therapeutics industries for the treatment of various metabolic diseases like type 2 diabetes and obesity. The structural–functional relationships of small molecule agonists and GLP-1R are yet to be understood. Therefore, an attempt was made on structurally known GLP-1R agonists (Compound 1, Compound 2, Compound A, Compound B, and (S)-8) to study their interaction with the extracellular domain of GLP-1R. In this study, we explored the dynamics, intrinsic stability, and binding mechanisms of these molecules through computational modeling, docking, molecular dynamics (MD) simulations and molecular mechanics Poisson–Boltzmann surface area (MM/PBSA) binding free energy estimation. Molecular docking study depicted that hydrophobic interaction (pi–pi stacking) plays a crucial role in maintaining the stability of the complex, which was also supported by intermolecular analysis from MD simulation study. Principal component analysis suggested that the terminal ends along with the turns/loops connecting adjacent helix and strands exhibit a comparatively higher movement of main chain atoms in most of the complexes. MM/PBSA binding free energy study revealed that non-polar solvation (van der Waals and electrostatic) energy subsidizes significantly to the total binding energy, and the polar solvation energy opposes the binding agonists to GLP-1R. Overall, we provide structural features information about GLP-1R complexes that would be conducive for the discovery of new GLP-1R agonists in the future for the treatment of various metabolic diseases.

Communicated by Ramaswamy H. Sarma     

Non-Boltzmann thermodynamic integration (NBTI) for macromolecular systems: relative free energy of binding of trypsin to benzamidine and benzylamine

The relative free energies of binding of trypsin to two amine inhibitors, benzamidine (BZD) and benzylamine (BZA), were calculated using non-Boltzmann thermodynamic integration (NBTI). Comparison of the simulations with the crystal structures of both complexes, trypsin-BZD and trypsin-BZA, shows that NBTI simulations better sample conformational space relative to thermodynamic integration (TI) simulations. The relative binding free energy calculated using NBTI was much closer to the experimentally determined value than that obtained using TI. The error in the TI simulation was found to be primarily due to incorrect sampling of BZA's conformation in the binding pocket. In contrast, NBTI produces a smooth mutation from BZD to BZA using a surrogate potential, resulting in a much closer agreement between the inhibitors' conformations and the omit electron density maps. This superior agreement between experiment and simulation, of both relative binding free energy differences and conformational sampling, demonstrates NBTI's usefulness for free energy calculations in macromolecular simulations.     

Probing the structural and energetic basis of kinesin-microtubule binding using computational alanine-scanning mutagenesis

Kinesin-microtubule (MT) binding plays a critical role in facilitating and regulating the motor function of kinesins. To obtain a detailed structural and energetic picture of kinesin-MT binding, we performed large-scale computational alanine-scanning mutagenesis based on long-time molecular dynamics (MD) simulations of the kinesin-MT complex in both ADP and ATP states. First, we built three all-atom kinesin-MT models: human conventional kinesin bound to ADP and mouse KIF1A bound to ADP and ATP. Then, we performed 30 ns MD simulations followed by kinesin-MT binding free energy calculations for both the wild type and mutants obtained after substitution of each charged residue of kinesin with alanine. We found that the kinesin-MT binding free energy is dominated by van der Waals interactions and further enhanced by electrostatic interactions. The calculated mutational changes in kinesin-MT binding free energy are in excellent agreement with results of an experimental alanine-scanning study with a root-mean-square error of ~0.32 kcal/mol [Woehlke, G., et al. (1997) Cell 90, 207-216]. We identified a set of important charged residues involved in the tuning of kinesin-MT binding, which are clustered on several secondary structural elements of kinesin (including well-studied loops L7, L8, L11, and L12, helices α4, α5, and α6, and less-explored loop L2). In particular, we found several key residues that make different contributions to kinesin-MT binding in ADP and ATP states. The mutations of these residues are predicted to fine-tune the motility of kinesin by modulating the conformational transition between the ADP state and the ATP state of kinesin.     

Binding free energy calculations between bovine β‐lactoglobulin and four fatty acids using the MMGBSA method

The bovine dairy protein β‐lactoglobulin (βlg) is a promiscuous protein that has the ability to bind several hydrophobic ligands. In this study, based on known experimental data, the dynamic interaction mechanism between bovine βlg and four fatty acids was investigated by a protocol combining molecular dynamics (MD) simulations and molecular mechanics generalized Born surface area (MMGBSA) binding free energy calculations. Energetic analyses revealed binding free energy trends that corroborated known experimental findings; larger ligand size corresponded to greater binding affinity. Finally, binding free energy decomposition provided detailed information about the key residues stabilizing the complex. © 2014 Wiley Periodicals, Inc. Biopolymers 101: 1010–1018, 2014.     

Absolute binding free energy calculations using molecular dynamics simulations with restraining potentials

The absolute (standard) binding free energy of eight FK506-related ligands to FKBP12 is calculated using free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent. A number of features are implemented to improve the accuracy and enhance the convergence of the calculations. First, the absolute binding free energy is decomposed into sequential steps during which the ligand-surrounding interactions as well as various biasing potentials restraining the translation, orientation, and conformation of the ligand are turned "on" and "off." Second, sampling of the ligand conformation is enforced by a restraining potential based on the root mean-square deviation relative to the bound state conformation. The effect of all the restraining potentials is rigorously unbiased, and it is shown explicitly that the final results are independent of all artificial restraints. Third, the repulsive and dispersive free energy contribution arising from the Lennard-Jones interactions of the ligand with its surrounding (protein and solvent) is calculated using the Weeks-Chandler-Andersen separation. This separation also improves convergence of the FEP/MD calculations. Fourth, to decrease the computational cost, only a small number of atoms in the vicinity of the binding site are simulated explicitly, while all the influence of the remaining atoms is incorporated implicitly using the generalized solvent boundary potential (GSBP) method. With GSBP, the size of the simulated FKBP12/ligand systems is significantly reduced, from approximately 25,000 to 2500. The computations are very efficient and the statistical error is small ( approximately 1 kcal/mol). The calculated binding free energies are generally in good agreement with available experimental data and previous calculations (within approximately 2 kcal/mol). The present results indicate that a strategy based on FEP/MD simulations of a reduced GSBP atomic model sampled with conformational, translational, and orientational restraining potentials can be computationally inexpensive and accurate.     

Molecular basis of the recognition of arachidonic acid by cytochrome P450 2E1 along major access tunnel

Cytochrome P450 2E1 is widely known for its ability to oxidize both low molecular weight xenobiotics and endogenous fatty acids (e.g., arachidonic acid (AA)). In this study, we investigated the structural features of the AA‐bound CYP2E1 complex utilizing molecular dynamics (MD) and found that the distinct binding modes for both AA and fatty acid analog are conserved. Moreover, multiple random acceleration MD simulations and steered MD simulations uncovered the most possible tunnel for fatty acids. The main attractions are derived from three key residues, His107, Ala108, and His109, whose side chains reorient to keep ligands bound via hydrogen bonds during the initial unbinding process. More importantly, based on the calculated binding free energy results, we hypothesize that the hydrogen bonds between the receptor and the ligand are the most important contributors involved in the binding affinity. Thus, it is inferred that the hydrogen bonds between these three residues and the ligand may help offer insights into the structural basis of the different ligand egress mechanisms for fatty acids and small weight compounds. Our investigation provides detailed atomistic insights into the structural features of human CYP2E1–fatty acid complex structures. Furthermore, the ligand‐binding characteristics obtained in the present study are helpful for both experimental and computational studies of CYPs and may allow future researchers to achieve desirable changes in enzymatic activities. © 2014 Wiley Periodicals, Inc. Biopolymers 103: 53–66, 2015.     

Viper venom hyaluronidase and its potential inhibitor analysis: a multipronged computational investigation

Viper venom hyaluronidase (VV-HYA) inhibitors have long been used as therapeutic agents for arresting the local and systemic effects caused during its envenomation. Henceforth, to understand its structural features and also to identify the best potential inhibitor against it the present computational study was undertaken. Structure-based homology modeling of VV-HYA followed by its docking and free energy-based ranking analysis of ligand, the MD simulations of the lead complex was also performed. The sequence analysis and homology modeling of VV-HYA revealed a distorted (β/α)₈ folding as in the case of hydrolases family of proteins. Molecular docking of the resultant 3D structure of VV-HYA with known inhibitors (compounds 1–25) revealed the importance of molecular recognition of hotspot residues (Tyr 75, Arg 288, and Trp 321) other than that of the active site residues. It also revealed that Trp 321 of VV-HYA is highly important for mediating π–π interactions with ligands. In addition, the molecular docking and comparative free energy binding analysis was investigated for the VV-HYA inhibitors (compounds 1–25). Both molecular docking and relative free energy binding analysis clearly confirmed the identification of sodium chromoglycate (compound 1) as the best potential inhibitor against VV-HYA. Molecular dynamics simulations additionally confirmed the stability of their binding interactions. Further, the information obtained from this work is believed to serve as an impetus for future rational designing of new novel VV-HYA inhibitors with improved activity and selectivity.     

Stabilization of the integrase‐DNA complex by Mg2+ ions and prediction of key residues for binding HIV‐1 integrase inhibitors

The HIV‐1 integrase is an attractive target for the therapeutics development against AIDS, as no host homologue of this protein has been identified. The integrase strand transfer inhibitors (INSTIs), including raltegravir, specifically target the second catalytic step of the integration process by binding to the DDE motif of the catalytic site and coordinating Mg²⁺ ions. Recent X‐ray crystallographic structures of the integrase/DNA complex from prototype foamy virus allowed to investigate the role of the different partners (integrase, DNA, Mg²⁺ ions, raltegravir) in the complex stability using molecular dynamics (MD) simulations. The presence of Mg²⁺ ions is found to be essential for the stability, whereas the simultaneous presence of raltegravir and Mg²⁺ ions has a destabilizing influence. A homology model of HIV‐1 integrase was built on the basis of the X‐ray crystallographic information, and protein marker residues for the ligand binding were detected by clustering the docking poses of known HIV‐1 integrase inhibitors on the model. Interestingly, we had already identified some of these residues to be involved in HIV‐1 resistance mutations and in the stabilization of the catalytic site during the MD simulations. Classification of protein conformations along MD simulations, as well as of ligand docking poses, was performed by using an original learning method, based on self‐organizing maps. This allows us to perform a more in‐depth investigation of the free‐energy basins populated by the complex in MD simulations on the one hand, and a straightforward classification of ligands according to their binding residues on the other hand. Proteins 2014; 82:466–478. © 2013 Wiley Periodicals, Inc.     

Thermodynamic basis of selectivity in guide-target-mismatched RNA interference

Silencing in RNAi is strongly affected by guide‐strand/target‐mRNA mismatches. Target nucleation is thought to occur at positions 2–8 of the guide (“seed region”); successful hybridization in this region is the primary determinant of target‐binding affinity and hence target cleavage. To define a molecular basis for the target sequence selectivity in RNAi, we studied all possible distinct single mismatches in seven positions of the seed region—a total of 21 substitutions. We report results from soft‐core thermodynamic integration simulations to determine changes in targeting binding‐free energies to Argonaute due to single mismatches in the guide strand, which arise during binding of an imperfectly matched target mRNA. In agreement with experiment, most mismatches impair target binding, consistent with a prominent role for binding affinity changes in RNAi sequence selectivity. Individual Argonaute residues located near the mismatched base pair are found to contribute significantly to binding affinity changes. We also use this methodology to analyze the mismatch‐dependent free energy changes for dissociation of a DNA?RNA hybrid from Argonaute, as a model for the escape of miRNAs from the silencing pathway. Several mismatched sequences of the miRNA have increased affinity to Argonaute, implying that some mismatches may reduce the probability for escape. Furthermore, calculations of base‐substitution‐dependent free energy changes for binding ssDNA reveal mild sequence sensitivity as expected for guide strand binding to Argonaute. Our findings give a thermodynamic basis for RNAi target sequence selectivity and suggest that miRNA mismatches may increase silencing effectiveness and thus could be evolutionarily advantageous. Proteins 2012; © 2011 Wiley Periodicals, Inc.     

Affinity enhancement of nanobody binding to EGFR: in silico site-directed mutagenesis and molecular dynamics simulation approaches

Epidermal growth factor receptor (EGFR), a transmembrane glycoprotein, is overexpressed in many cancers such as head-neck, breast, prostate, and skin cancers for this reason it is a good target in cancer therapy and diagnosis. In nanobody-based cancer diagnosis and treatment, nanobodies with high affinity toward receptor (e.g. EGFR) results in effective treatment or diagnosis of cancer. In this regard, the main aim of this study is to develop a method based on molecular dynamic (MD) simulations for designing of 7D12 based nanobody with high affinity compared with wild-type nanobody. By surveying electrostatic and desolvation interactions between different residues of 7D12 and EGFR, the critical residues of 7D12 that play the main role in the binding of 7D12 to EGFR were elucidated and based on these residues, five logical variants were designed. Following the 50 ns MD simulations, pull and umbrella sampling simulation were performed for 7D12 and all its variants in complex with EGFR. Binding free energy of 7D12 (and all its variants) with EGFR was obtained by weighted histogram analysis method. According to binding free energy results, GLY101 to GLU mutation showed the highest binding affinity but this variant is unstable after 50 ns MD simulations. ALA100 to GLU mutation shows suitable binding enhancement with acceptable structural stability. Suitable position and orientation of GLU in residue 100 of 7D12 against related amino acids of EGFR formed some extra hydrogen and electrostatic interactions which resulted in binding enhancement.     

Improving the accuracy of protein stability predictions with multistate design using a variety of backbone ensembles

Multistate computational protein design (MSD) with backbone ensembles approximating conformational flexibility can predict higher quality sequences than single‐state design with a single fixed backbone. However, it is currently unclear what characteristics of backbone ensembles are required for the accurate prediction of protein sequence stability. In this study, we aimed to improve the accuracy of protein stability predictions made with MSD by using a variety of backbone ensembles to recapitulate the experimentally measured stability of 85 Streptococcal protein G domain β1 sequences. Ensembles tested here include an NMR ensemble as well as those generated by molecular dynamics (MD) simulations, by Backrub motions, and by PertMin, a new method that we developed involving the perturbation of atomic coordinates followed by energy minimization. MSD with the PertMin ensembles resulted in the most accurate predictions by providing the highest number of stable sequences in the top 25, and by correctly binning sequences as stable or unstable with the highest success rate (≈90%) and the lowest number of false positives. The performance of PertMin ensembles is due to the fact that their members closely resemble the input crystal structure and have low potential energy. Conversely, the NMR ensemble as well as those generated by MD simulations at 500 or 1000 K reduced prediction accuracy due to their low structural similarity to the crystal structure. The ensembles tested herein thus represent on‐ or off‐target models of the native protein fold and could be used in future studies to design for desired properties other than stability. Proteins 2014; 82:771–784. © 2013 Wiley Periodicals, Inc.