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.     

Exploration of the chlorpyrifos escape pathway from acylpeptide hydrolases using steered molecular dynamics simulations

Acylpeptide hydrolases (APH) catalyze the removal of an N-acylated amino acid from blocked peptides. APH is significantly more sensitive than acetylcholinesterase, a target of Alzheimer’s disease, to inhibition by organophosphorus (OP) compounds. Thus, OP compounds can be used as a tool to probe the physiological functions of APH. Here, we report the results of a computational study of molecular dynamics simulations of APH bound to the OP compounds and an exploration of the chlorpyrifos escape pathway using steered molecular dynamics (SMD) simulations. In addition, we apply SMD simulations to identify potential escape routes of chlorpyrifos from hydrolase hydrophobic cavities in the APH-inhibitor complex. Two previously proposed APH pathways were reliably identified by CAVER 3.0, with the estimated relative importance of P1 > P2 for its size. We identify the major pathway, P2, using SMD simulations, and Arg526, Glu88, Gly86, and Asn65 are identified as important residues for the ligand leaving via P2. These results may help in the design of APH-targeting drugs with improved efficacy, as well as in understanding APH selectivity of the inhibitor binding in the prolyl oligopeptidase family.     

The mechanism for the complexation and dissociation between siRNA and PMAL: a molecular dynamics simulation study based on a coarse-grained model

Abstract

The capacity of silencing genes makes small interfering RNA (siRNA) becomes potential candidates for curing many fatal diseases. Due to the low stability and delivery efficiency of siRNA, the design of amphiphilic carrier for siRNA delivery is vital for the practical gene therapy. In the present work, we explored how the complexation and dissociation of siRNA with poly (maleic anhydride-alt-1-decene) substituted with 3-(dimethylamino) propylamine (PMAL), which is a recent synthesised amphiphilic polymer and can be used in delivery of siRNAs and proteins, using traditional molecular dynamics simulations, together with steered molecular dynamics simulations. It was shown that the complexation of siRNA with PMALs can spontaneously occur, no matter what unit number of PMAL is. PMALs of different unit numbers form micelle-like structures and interact with siRNA surface. With the increase of unit number, PMAL becomes more flexible and interacts with siRNA from attachment to entanglement. The dissociation of PMAL from siRNA is an energy-consuming process. The free energy difference increases with the unit number of PMAL. The free energy for dissociation involves both the stretch of PMAL and the separation of PMAL from siRNA. Therefore, an optimal unit number of PMAL is critical for the delivery efficiency of siRNA when PMAL is used as carrier. In present work, when the radius of gyration of PMAL approaches to that of siRNA, PMAL gives a favoured both complexation and dissociation between siRNA and PMAL. Finally, we propose the mechanism of complexation and dissociation of PMAL with siRNA. The above simulation established a molecular insight of the interaction between siRNA and PMAL and was helpful for the design and applications of new PMAL-based polymers as siRNA delivery carriers.     

Ligand unbinding pathways from the vitamin D receptor studied by molecular dynamics simulations

Molecular dynamics simulation techniques have been used to study the unbinding pathways of 1α,25-dihydroxyvitamin D₃ from the ligand-binding pocket of the vitamin D receptor (VDR). The pathways observed in a large number of relatively short (<200 ps) random acceleration molecular dynamics (RAMD) trajectories were found to be in fair agreement, both in terms of pathway locations and deduced relative preferences, compared to targeted molecular dynamics (TMD) and streered molecular dynamics simulations (SMD). However, the high-velocity ligand expulsions of RAMD tend to favor straight expulsion trajectories and the observed relative frequencies of different pathways were biased towards the probability of entering a particular exit channel. Simulations indicated that for VDR the unbinding pathway between the H1–H2 loop and the β-sheet between H5 and H6 is more favorable than the pathway located between the H1–H2 loop and H3. The latter pathway has been suggested to be the most likely unbinding path for thyroid hormone receptors (TRs) and a likely path for retinoic acid receptor. Ligand entry/exit through these two pathways would not require displacement of H12 from its agonistic position. Differences in the packing of the H1, H2, H3 and β-sheet region explain the changed relative preference of the two unbinding pathways in VDR and TRs. Based on the crystal structures of the ligand binding domains of class 2 nuclear receptors, whose members are VDR and TRs, this receptor class can be divided in two groups according to the packing of the H1, H2, H3 and β-sheet region. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Insight into the interactive residues between two domains of human somatic Angiotensin-converting enzyme and Angiotensin II by MM-PBSA calculation and steered molecular dynamics simulation

Angiotensin-converting enzyme (ACE), a membrane-bound zinc metallopeptidase, catalyzes the formation of Angiotensin-II (AngII) and the deactivation of bradykinin in the renin–angiotensin-aldosterone and kallikrein–kinin systems. As a hydrolysis product of ACE, AngII is regarded as an inhibitor and displays stronger competitive inhibition in the C-domain than the N-domain of ACE. However, the AngII binding differences between the two domains and the mechanisms behind AngII dissociation from the C-domain are rarely explored. In this work, molecular docking, Molecular Mechanics/Poisson–Boltzmann Surface Area calculation, and steered molecular dynamics (SMD) are applied to explore the structures and interactions in the binding or unbinding of AngII with the two domains of human somatic ACE. Calculated free energy values suggest that the C-domain–AngII complex is more stable than the N-domain–AngII complex, consistent with available experimental data. SMD simulation results imply that electrostatic interaction is dominant in the dissociation of AngII from the C-domain. Moreover, Gln106, Asp121, Glu123, and Tyr213 may be the key residues in the unbinding pathway of AngII. The simulation results in our work provide insights into the interactions between the two domains of ACE and its natural peptide inhibitor AngII at a molecular level. Moreover, the results provide theoretical clues for the design of new inhibitors.     

Critical role of helix 12 of the vitamin D(3) receptor for the partial agonism of carboxylic ester antagonists.

The carboxy-terminal alpha-helix of a nuclear receptor ligand-binding domain (LBD), helix 12, contains a critical, ligand-modulated interface for the interaction with coactivator proteins. In this study, using the example of the vitamin D receptor (VDR) and the partial antagonist ZK159222, the role of helix 12 (residues 417-427) for both antagonistic and agonistic receptor actions was investigated. Amino acid residue G423 was demonstrated to be critical for partial agonism of ZK159222, but not for the activity of the natural VDR agonist, 1alpha,25-dihydroxyvitamin D(3) (1alpha,25(OH)(2)D(3)). The amount of partial agonism of ZK159222 increased when helix 12 was truncated by the last four amino acid residues (Delta424-27) and augmented even more, when in addition helix 12 of VDR's dimerization partner, retinoid X receptor (RXR), was truncated. In contrast, the low agonism of a structural derivative of ZK159222, ZK168281, was not affected comparably, whereas other close structural relatives of ZK159222 even demonstrated the same agonistic activity as that of 1alpha,25(OH)(2)D(3). The amount of agonism of ZK159222 and ZK168281 at different variations of helix 12 correlated well with VDR's ability to complex with coactivator proteins and inversely correlated with the strength of the compound's antagonistic action on 1alpha,25(OH)(2)D(3) signalling. Molecular dynamics simulations of the LBD complexed with the two antagonists could explain their different action by demonstrating a more drastic displacement of helix 12 through ZK168281 than through ZK159222. Moreover, the modelling could indicate a kink of helix 12 at amino acid residue G423, which provides the last four amino acid residues of helix 12 with a modulatory role for the partial agonism of some VDR antagonists, such as ZK159222. In conclusion, partial agonism of a VDR antagonist is lower the more it disturbs helix 12 in taking the optimal position for coactivator interaction.     

Insight into the process of product expulsion in cellobiohydrolase Cel6A from Trichoderma reesei by computational modeling

Glycoside hydrolase cellulase family 6 from Trichoderma reesei (TrCel6A) is an important cellobiohydrolase to hydrolyze cellooligosaccharide into cellobiose. The knowledge of enzymatic mechanisms is critical for improving the conversion efficiency of cellulose into ethanol or other chemicals. However, the process of product expulsion, a key component of enzymatic depolymerization, from TrCel6A has not yet been described in detail. Here, conventional molecular dynamics and steered molecular dynamics (SMD) were applied to study product expulsion from TrCel6A. Tyr103 may be a crucial residue in product expulsion given that it exhibits two different posthydrolytic conformations. In one conformation, Tyr103 rotates to open the ?3 subsite. However, Tyr103 does not rotate in the other conformation. Three different routes for product expulsion were proposed on the basis of the two different conformations. The total energy barriers of the three routes were calculated through SMD simulations. The total energy barrier of product expulsion through Route 1, in which Tyr103 does not rotate, was 22.2 kcal·mol^?1. The total energy barriers of product expulsion through Routes 2 and 3, in which Tyr103 rotates to open the ?3 subsite, were 10.3 and 14.4 kcal·mol^?1, respectively. Therefore, Routes 2 and 3 have lower energy barriers than Route 1, and Route 2 is the thermodynamically optimal route for product expulsion. Consequently, the rotation of Tyr103 may be crucial for product release from TrCel6A. Results of this work have potential applications in cellulase engineering.     

p53 DNA结合域与抗癌多肽的拉伸分子动力学研究

研究发现,取自蓝铜蛋白azurin的一段多肽p28能够进入癌细胞,结合到肿瘤抑制因子p53的DNA结合域上,进而增加p53的抗癌能力.本工作中,通过拉伸分子动力学方法,在原子尺度上研究了p28-p53 DBD复合物的解离过程.分析结果显示复合物的解离过程遵循着一定的分离顺序.对解离力的分析以及对沿着解离路径的不可逆做功的计算,使我们能够从复合物的能量地貌中提取有用的信息,而这些信息也决定了复合物的解离过程.     

Drug repurposing for chronic myeloid leukemia: in silico and in vitro investigation of DrugBank database for allosteric Bcr-Abl inhibitors

Chronic myeloid leukemia (CML) is caused by chromosomal rearrangement resulting in the expression of Bcr-Abl fusion protein with deregulated Abl tyrosine kinase activity. Approved drugs – imatinib, dasatinib, nilotinib, and ponatinib – target the ATP-binding site of Abl kinase. Even though these drugs are initially effective, long-term usefulness is limited by the development of resistance. To overcome this problem, targeting the allosteric site of Abl kinase, which is remote from the ATP-binding site is found to be a useful strategy. In this study, structure-based and ligand-based virtual screening methods were applied to narrow down possible drugs (from DrugBank database) that could target the allosteric site of Abl kinase. Detailed investigations of the selected drugs in the allosteric site of Abl kinase, using molecular dynamics and steered molecular dynamics simulation shows that gefitinib, an EGFR inhibitor approved for the treatment of lung cancer, could bind effectively to the allosteric site of Bcr-Abl. More interestingly, gefitinib was found to enhance the ability of imatinib to bind at the ATP-binding site of Bcr-Abl kinase. Based on the in silico findings, gefitinib was tested in combination with imatinib in K562 CML cell line using MTT cell proliferation assay and found to have a synergistic antiproliferative activity. Further detailed mechanistic study could help to unravel the full potential of imatinib – gefitinib combination for the treatment of CML.     

Donor-strand exchange in chaperone-assisted pilus assembly revealed in atomic detail by molecular dynamics

Adhesive multi-subunit fibres are assembled on the surface of many pathogenic bacteria via the chaperone-usher pathway. In the periplasm, a chaperone donates a β-strand to a pilus subunit to complement its incomplete immunoglobulin-like fold. At the outer membrane, this is replaced with a β-strand formed from the N-terminal extension (Nte) of an incoming pilus subunit by a donor-strand exchange (DSE) mechanism. This reaction has previously been shown to proceed via a concerted mechanism, in which the Nte interacts with the chaperone:subunit complex before the chaperone has been displaced, forming a ternary intermediate. Thereafter, the pilus and chaperone β-strands have been postulated to undergo a strand swap by a ‘zip-in-zip-out’ mechanism, whereby the chaperone strand zips out, residue by residue, as the Nte simultaneously zips in, although direct experimental evidence for a zippering mechanism is still lacking. Here, molecular dynamics simulations have been used to probe the DSE mechanism during formation of the Saf pilus from Salmonella enterica at the atomic level, allowing the direct investigation of the zip-in-zip-out hypothesis. The simulations provide an explanation of how the incoming Nte is able to dock and initiate DSE due to inherent dynamic fluctuations within the chaperone:subunit complex. In the simulations, the chaperone donor strand was seen to unbind from the pilus subunit, residue by residue, in direct support of the zip-in-zip-out hypothesis. In addition, an interaction of a residue towards the N-terminus of the Nte with a specific binding pocket (P*) on the adjacent pilus subunit was seen to stabilise the DSE product against unbinding, which also proceeded in the simulations by a zippering mechanism. Together, the study provides an in-depth picture of DSE, including the first atomistic insights into the molecular events occurring during the zip-in-zip-out mechanism.     

Exploring stereochemical specificity of phosphotriesterase by MM-PBSA and MM-GBSA calculation and steered molecular dynamics simulation

Wild-type phosphotriesterase (PTE) prefers the S_P-enantiomers over the corresponding R_P-enantiomers by factors ranging from 10 to 90. To satisfy the binding modes of the PTE of S_P- and R_P-enantiomers, all-atom molecular dynamics simulations were carried out on two paraoxon S_P and R_P derivatives, namely, S_p-1 and R_p-1. Molecular mechanics Poisson–Boltzmann surface area and molecular mechanics generalized Born surface area (MM-PBSA and MM-GBSA) calculations indicated that His230 in S_p-1-PTE had a closer interaction with the substrate than that in R_p-1-PTE and that such interaction increased the catalytic efficiency of PTE for S_p-1. The steered molecular dynamics simulation indicated that, compared with S_p-1, R_p-1 in the unbinding (binding) may hinder some residue displacement, thus requiring more effort to escape the binding pocket of PTE. In addition, Trp131, Phe306, and Tyr309 are deemed important residues for the S_p-1 unbinding pathway via PTE, whereas Tyr309 alone is considered an important residue for the R_p-1 unbinding pathway. These results demonstrate the possibility of dramatically altering the stereoselectivity and overall reactivity of the native enzyme toward chiral substrates by modifying specific residues located within the active site of PTE.     

Investigations of Takeout proteins’ ligand binding and release mechanism using molecular dynamics simulation

Takeout (To) proteins exist in a diverse range of insect species. They are involved in many important processes of insect physiology and behaviors. As the ligand carriers, To proteins can transport the small molecule to the target tissues. However, ligand release mechanism of To proteins is unclear so far. In this contribution, the process and pathway of the ligand binding and release are revealed by conventional molecular dynamics simulation, steered molecular dynamics simulation and umbrella sampling methods. Our results show that the α4-side of the protein is the unique gate for the ligand binding and release. The structural analysis confirms that the internal cavity of the protein has high rigidity, which is in accordance with the recent experimental results. By using the potential of mean force calculations in combination with residue cross correlation calculation, we concluded that the binding between the ligand and To proteins is a process of conformational selection. Furthermore, the conformational changes of To proteins and the hydrophobic interactions both are the key factors for ligand binding and release.     

Pentameric ligand-gated ion channels: insights from computation

Pentameric ligand-gated ion channels (pLGICs) conduct upon the binding of an agonist and are fundamental to neurotransmission. New insights into the complex mechanisms underlying pLGIC gating, ion selectivity and modulation have recently been gained via a series of crystal structures in prokaryotes and Caenorhabditis elegans, as well as computational studies relying on these structures. Here, we review contributions from a variety of computational approaches, including normal-mode analysis, automated docking and fully atomistic molecular dynamics simulation. Examples from our own research, particularly concerning interactions with general anaesthetics and lipids, are used to illustrate predictive results complementary to crystallographic studies.