A two-state homology model of the hERG K+ channel: application to ligand binding

Homology models based on available K+ channel structures have been used to construct a multiple state representation of the hERG cardiac K+ channel. These states are used to capture the flexibility of the channel. We show that this flexibility is essential in order to correctly model the binding affinity of a set of diverse ligands. Using this multiple state approach, a binding affinity model was constructed for set of known hERG channel binders. The predicted pIC50s are in good agreement with experiment (RMSD: 0.56 kcal/mol). In addition, these calculations provide structures for the bound ligands that are consistent with published mutation studies. These computed ligand bound complex structures can be used to guide synthesis of analogs with reduced hERG liability.     

Insight into the molecular interaction between the cyclic nucleotide-binding homology domain and the eag domain of the hERG channel

The gating of the hERG channel is regulated by its eag domain through molecular interaction with either the cyclic nucleotide-binding homology domain (CNBHD) or the linker between transmembrane segments 4 and 5. Our NMR study on the purified CNBHD demonstrated that it contains nine β-strands and does not bind cAMP. We show that the eag domain binds to the CBND through an interface containing several disease-associated mutations. The N-terminal cap domain and R56 in the eag domain are important for the interaction with the CNBHD. Residues from the CNBHD that were affected by the interaction with the eag domain were also identified. A R56Q mutation does not cause major structural changes in the eag domain and showed reduced interaction with the CNBHD.     

Development of a potent 2-oxoamide inhibitor of secreted phospholipase A2 guided by molecular docking calculations and molecular dynamics simulations

Inhibition of group IIA secreted phospholipase A₂ (GIIA sPLA₂) has been an important objective for medicinal chemists. We have previously shown that inhibitors incorporating the 2-oxoamide functionality may inhibit human and mouse GIIA sPLA₂s. Herein, the development of new potent inhibitors by molecular docking calculations using the structure of the known inhibitor 7 as scaffold, are described. Synthesis and biological evaluation of the new compounds revealed that the long chain 2-oxoamide based on (S)-valine GK241 led to improved activity (IC₅₀ = 143 nM and 68 nM against human and mouse GIIA sPLA₂, respectively). In addition, molecular dynamics simulations were employed to shed light on GK241 potent and selective inhibitory activity.     

A specific two-pore domain potassium channel blocker defines the structure of the TASK-1 open pore

Two-pore domain potassium (K(2P)) channels play a key role in setting the membrane potential of excitable cells. Despite their role as putative targets for drugs and general anesthetics, little is known about the structure and the drug binding site of K(2P) channels. We describe A1899 as a potent and highly selective blocker of the K(2P) channel TASK-1. As A1899 acts as an open-channel blocker and binds to residues forming the wall of the central cavity, the drug was used to further our understanding of the channel pore. Using alanine mutagenesis screens, we have identified residues in both pore loops, the M2 and M4 segments, and the halothane response element to form the drug binding site of TASK-1. Our experimental data were used to validate a K(2P) open-pore homology model of TASK-1, providing structural insights for future rational design of drugs targeting K(2P) channels.     

Insights into structure and function of SHIP2-SH2: homology modeling, docking, and molecular dynamics study

SRC homology 2 (SH2)-containing inositol 5′-phosphatase protein (SHIP2) is a potential target for type 2 diabetes. Its ability to dephosphorylate the lipid messenger phosphatidylinositol 3,4,5-trisphosphate [PtdIns(3,4,5)P3], important for insulin signaling, makes it an important target against type 2 diabetes. The insulin-induced SHIP2 interaction with Shc is very important for the membrane localization and functioning of SHIP2. There is a bidentate relationship between the two proteins where two domains each from SHIP2 and Shc are involved in mutual binding. However in the present study, the SHIP2-SH2 domain binding with the phosphorylated tyrosine 317 on the collagen-homology (CH) domain of Shc, has been studied due to the indispensability of this interaction in SHIP2 localization. In the absence of the crystal structure of SHIP2-SH2, its structural model was developed followed by tracking its molecular interactions with Shc through molecular docking and dynamics studies. This study revealed much about the structural interactions between the SHIP2-SH2 and Shc-CH. Finally, docking study of a nonpeptide inhibitor into the SHIP2-SH2 domain further confirmed the structural interactions involved in ligand binding and also proposed the inhibitor as a major starting point against SHIP2-SH2 inhibition. The insights gained from the current study should prove useful in the design of more potent inhibitors against type 2 diabetes.     

On the putative binding site of RFamide-family neuropeptides from the western Atlantic clam Sunray Venus and cephalopods on acid-sensing ion channels. An automated docking and molecular-dynamics study with hASIC1a homology model

Investigated here are interactions of C-terminal amidated peptides with the hASIC1a acid-sensing ion channel. The peptides comprise endogenous FMRFa, present in the western Atlantic clam Sunray Venus, and FIRFa, present in cephalopods, as well as non-endogenous ones for comparison. The interaction is investigated by automated docking. The resulting key hASIC1a-FMRFa complex, set in a lipidic POPC (=1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) membrane surrounded by H(2) O and Na(+) -neutralized, was also investigated by molecular dynamics. It was observed that all investigated peptides become encapsulated into the ion channel, on one side by the thumb and finger of a subunit, and, on the opposite side, by the knuckle and β-ball of a second subunit. The third subunit is not involved. This is much the same binding site that was disclosed previously by both a similar computational approach, and electrophysiological and binding experiments for the hASIC1a ion channel-blocker tarantula toxin PCTX1. This paves the way to a better understanding of the role of these peptides in invertebrates.     

Combined docking and molecular dynamics simulations to enlighten the capacity of Pseudomonas cepacia and Candida antarctica lipases to catalyze quercetin acetylation

A combined docking and molecular dynamics protocol was applied to investigate quercetin binding modes within the catalytic cavity of Candida antarctica lipase B (CALB) and Pseudomonas cepacia lipase (PCL), aiming to explain the difference of specificity of these enzymes in acetylation reaction. For both lipases, docking of quercetin yielded two families of conformers with either the quercetin A or B-ring pointing towards the catalytic residues. Molecular dynamics (MD) calculations were subsequently performed on several complexes of each family. MD trajectories were analyzed focusing on the orientation of the acyl donor bound to the catalytic serine towards the oxyanion hole residues and the proximity of quercetin hydroxyl groups to the catalytic residues. Results showed that with CALB, the acetate was not correctly positioned within the oxyanion hole whatever the orientation of quercetin, suggesting that no product could be obtained. With PCL, the acetate remained within the oxyanion hole during all MD trajectories. Depending on quercetin orientation, either the 7-OH group or the 3, 5, 3′, 4′-OH groups came alternatively near the catalytic residues, suggesting that all of them could be acylated. The capacity of models to explain the regioselectivity of the reaction was discussed. Key residues and interactions involved in quercetin binding modes were identified and related to the reaction feasibility.     

Electrostatics studies and molecular dynamics simulations of a homology model of the Shaker K+ channel pore

A homology model of the pore domain of the Shaker K+ channel has been constructed using a bacterial K+ channel, KcsA, as a template structure. The model is in agreement with mutagenesis and sequence variability data. A number of structural features are conserved between the two channels, including a ring of tryptophan sidechains on the outer surface of the pore domain at the extracellular end of the helix bundle, and rings of acidic sidechains close to the extracellular mouth of the channel. One of these rings, that formed by four Asp447 sidechains at the mouth of the Shaker pore, is shown by pK(A) calculations to be incompletely ionized at neutral pH. The potential energy profile for a K+ ion moved along the central axis of the Shaker pore domain model selectivity filter reveals a shallow well, the depth of which is modulated by the ionization state of the Asp447 ring. This is more consistent with the high cation flux exhibited by the channel in its conductance value of 19 pS.     

Exploring the binding mode of donepezil with calf thymus DNA using spectroscopic and molecular docking methods

Donepezil (DNP) is one of approved drugs to treat Alzheimer's disease (AD). However, the potential effect of DNP on DNA is still unclear. Therefore, the interaction of DNP with calf thymus DNA (DNA) was studied in vitro using spectroscopic and molecular docking methods. Steady‐state and transient fluorescence experiments showed that there was a clear binding interaction between DNP and DNA, resulting from DNP fluorescence being quenched using DNA. DNP and DNA have one binding site between them, and the binding constant (K_b) was 0.78 × 10⁴ L·mol^?1 at 298 K. In this binding process, hydrophobic force was the main interaction force, because enthalpy change (ΔH) and entropy change (ΔS) of DNP–DNA were 67.92 kJ·mol^?1 and 302.96 J·mol^?1·K^?1, respectively. DNP bound to DNA in a groove‐binding mode, which was verified using a competition displacement study and other typical spectroscopic methods. Fourier transform infrared (FTIR) spectrum results showed that DNP interacted with guanine (G) and cytosine (C) bases of DNA. The molecular docking results further supported the results of spectroscopic experiments, and suggested that both Pi‐Sigma force and Pi‐Alkyl force were the major hydrophobic force functioning between DNP and DNA.     

Prediction of the binding modes between BB-83698 and peptide deformylase from Bacillus stearothermophilus by docking and molecular dynamics simulation

BB-83698 is a first potent inhibitor of peptide deformylase in this novel class to enter clinical trials. In this study, automated docking, molecular dynamics simulation and binding free energy calculations with the linear interaction energy (LIE) method are first applied to investigate the binding of BB-83698 to the peptide deformylase from Bacillus stearothermophilus. The lowest docking energy structure from each cluster is selected as different representative binding modes. Compared with the experimental data, the results show that the binding of BB-83698 in Mode 1 is the most stable, with a binding free energy of -41.35 kJ/mol. The average structure of the Mode 1 complex suggests that inhibitor interacts with Ile59 and Gly109 by hydrogen bond interaction and with Pro47, Pro57, Ile59 and Leu146 by hydrophobic interaction are essential for the activity of BB-83698. Mode 2 represents a new binding mode. Additionally, if the hydrophilic group is introduced to the benzo-[1,3]-dioxole ring, the binding affinity of BB-83698 to the peptide deformylase from B. stearothermophilus will be greatly improved.     

Ion binding properties and structure stability of the NaK channel

Ion distribution in the selectivity filter and ion-water and ion-protein interactions of NaK channel are systematically investigated by all-atom molecular dynamics simulations, with the tetramer channel protein being embedded in a solvated phospholipid bilayer. Analysis of the simulation results indicates that K⁺ ions prefer to bind within the sites formed by two adjacent planes of oxygen atoms from the selectivity filter, while Na⁺ ions are inclined to bind to a single plane of four oxygen atoms. At the same time, both K⁺ and Na⁺ ions can diffuse in the vestibule, accompanying with movements of the water molecules confined in a complex formed by the vestibule together with four small grottos connecting to it. As a result, K⁺ ions show a wide range of coordination numbers (6-8), while Na⁺ ions display a constant coordination number of ∼ 6 in the selectivity filter, which may result in the loss of selectivity of NaK. It is also found that a Ca²⁺ can bind at the extracellular site as reported in the crystal structure in a partially hydrated state, or at a higher site in a full hydration state. Furthermore, the carbonyl group of Asp66 can reorient to point towards the center pore when an ion exists in the vestibule, while that of Gly65 always aligns tangentially to the channel axis, as in the crystallographic structures.     

Investigating the disorder–order transition of calmodulin binding domain upon binding calmodulin using molecular dynamics simulation

We have studied the conformational transition of the calmodulin binding domains (CBD) in several calmodulin‐binding kinases, in which CBD changes from the disordered state to the ordered state when binding with calmodulin (CaM). Targeted molecular dynamics simulation was used to investigate the binding process of CaM and CBD of CaM‐dependent kinase I (CaMKI–CBD). The results show that CaMKI–CBD began to form an α‐helix and the interaction free energy between CaM and CaMKI–CBD increased once CaM fully encompassed CaMKI–CBD. Two series of CaM/CBD complex systems, including the complexes of CaM with the initially disordered and the final ordered CBD, were constructed to study the interaction using molecular dynamics simulations. Our analyses suggest that the VDW interaction plays a dominant role in CaM/CBD binding and is a key factor in the disorder–order transition of CBD. Additionally, the entropy effect is not in favor of the formation of the CaM/CBD complex, which is consistent with the experimental evidence. Based on the results, it appears that the CBD conformational change from a non‐compact extended structure to compact α‐helix is critical in gaining a favorable VDW interaction and interaction free energy. Copyright © 2009 John Wiley & Sons, Ltd.     

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.     

Interaction between K+ channel gate modifier hanatoxin and lipid bilayer membranes analyzed by molecular dynamics simulation

Hanatoxin (HaTx) is an ellipsoidal-shaped peptide that binds to the voltage sensor of voltage-dependent channels. Of physicochemical interest, HaTx has a “ring” of charged residues around its periphery and a hydrophobic protrusion. It has previously been postulated that HaTx binds to and functions on the surface of membranes, but a recent fluorescent-quenching study has implied a fairly deep positioning of HaTx in the lipid bilayer membrane. We carried out numerous molecular dynamic simulations of HaTx1, a well-studied variant of HaTx, in fully hydrated phospholipid bilayers. The system reproduced the surface-binding mode of HaTx1, in which HaTx1 resided in the extracellular side (outer) of the water/membrane interface with the hydrophobic patch of HaTx1 facing the membrane interior. On the other hand, analyses with various parameter settings suggested that the surface-binding mode was unstable because of the substantial attractive electrostatic force between HaTx1 and the lipid head groups of the inner (opposite) leaflet. Compared with this electrostatic force, the energetic cost for membrane deformation involving meniscus formation appeared to be small. In an attempt to interpret the quenching data, we consider the possibility of dimpling (meniscus formation) that brings HaTx1 inward (only ~0.7–0.8 nm above the bilayer center), while accounting for the flexibility of both leaflets of the membrane and the long-range interaction between positively charged residues of the membrane-bound peptide and the polar head groups of the opposite leaflet of the membrane. It is suggested that molecular dynamics simulations taking into account the flexibility of the membrane surface is potentially useful in interpreting the fluorescence-quenching data.Electronic Supplementary Material Supplementary material is available for this article at and is accessible for authorized users.     

From the Sequence to the Conformation of the Unabridged Transmembrane Domains TM1 and TM2 of the cASIC1a Ion Channel – A Parallel Tempering Approach

This work was devised to unravel, along replica‐exchange molecular‐dynamics (REMD) simulations, the conformation in solution of the TM1 and TM2 transmembrane domains of the homotrimeric cASIC1a ion channel. This includes the head of TM1 and tail of TM2 that had previously defied X‐ray diffraction analysis in the crystal. The structure of the open‐channel complex of cASIC1a with psalmotoxin 1 (PcTx1) was chosen here as a basis, although, to make the simulations affordable, the procedure was limited to the missing portions, including a few adjacent α‐helical turns. The latter were held fixed during the simulations. Reassembling the whole subunit, by superimposition of the fixed portions, resulted in diving of both TM1 and TM2 as continuous α‐helices into the cytoplasm. At completion of this work, it appeared, from similar X‐ray diffraction studies, that TM2 for both the complex of cASIC1a with the coral snake MitTx toxin, and the isolated desensitized ion channel, is discontinuous, with the triad G443‐A444‐S445 taking an extended, belt‐like conformation. In this way, a filter ring against hydrated ions is formed by G443 in the trimer. Our REMD examination of this complex revealed a strong resistance by G443, and only that residue, to take dihedral‐angle values compatible with an α‐helical conformation. This suggests that the flexibility of glycine alone does not explain formation of the extended, belt‐like conformation of the triad G443‐A444‐S445. This also requires cooperation in the trimer.     

Non-Ewald methods: theory and applications to molecular systems

Several non-Ewald methods for calculating electrostatic interactions have recently been developed, such as the Wolf method, the reaction field method, the pre-averaging method, and the zero-dipole summation method, for molecular dynamics simulations of various physical systems, including biomolecular systems. We review the theories of these approaches and their potential applications to molecular simulations, and discuss their relationships.     

Exploring the energy landscape of protein folding using replica-exchange and conventional molecular dynamics simulations

Two independent replica-exchange molecular dynamics (REMD) simulations with an explicit water model were performed of the Trp-cage mini-protein. In the first REMD simulation, the replicas started from the native conformation, while in the second they started from a nonnative conformation. Initially, the first simulation yielded results qualitatively similar to those of two previously published REMD simulations: the protein appeared to be over-stabilized, with the predicted melting temperature 50-150K higher than the experimental value of 315K. However, as the first REMD simulation progressed, the protein unfolded at all temperatures. In our second REMD simulation, which starts from a nonnative conformation, there was no evidence of significant folding. Transitions from the unfolded to the folded state did not occur on the timescale of these simulations, despite the expected improvement in sampling of REMD over conventional molecular dynamics (MD) simulations. The combined 1.42 micros of simulation time was insufficient for REMD simulations with different starting structures to converge. Conventional MD simulations at a range of temperatures were also performed. In contrast to REMD, the conventional MD simulations provide an estimate of Tm in good agreement with experiment. Furthermore, the conventional MD is a fraction of the cost of REMD and continuous, realistic pathways of the unfolding process at atomic resolution are obtained.     

Molecular modeling studies of atorvastatin analogues as HMGR inhibitors using 3D-QSAR,molecular docking and molecular dynamics simulations

3-Hydroxy-3-methylglutaryl coenzyme-A reductase (HMGR) is generally regarded as targets for the treatment of hypercholesterolemia. HMGR inhibitors (more commonly known as statins) are discovered as plasma cholesterol lowering molecules. In this work, 120 atorvastatin analogues were studied using a combination of molecular modeling techniques including three-dimensional quantitative structure–activity relationship (3D-QSAR), molecular docking and molecular dynamics (MD) simulation. The results show that the best CoMFA (comparative molecular field analysis) model has q² = 0.558 and r² = 0.977, and the best CoMSIA (comparative molecular similarity indices analysis) model has q² = 0.582 and r² = 0.919. Molecular docking and MD simulation explored the binding relationship of the ligand and the receptor protein. The calculation results indicated that the hydrophobic and electrostatic fields play key roles in QSAR model. After MD simulation, we found four vital residues (Lys735, Arg590, Asp690 and Asn686) and three hydrophobic regions in HMGR binding site. The calculation results show that atorvastatin analogues obtained by introduction of F atoms or gem-difluoro groups could obviously improve the inhibitory activity. The new HMGR inhibitor analogues design in this Letter had been submitted which is being currently synthesized by our laboratories.     

Determination on the binding of chlortetracycline to bovine serum albumin using spectroscopic methods

In this work, the interaction of chlortetracycline with bovine serum albumin (BSA) was investigated by fluorescence spectroscopy, circular dichroism (CD) spectroscopy, and molecular docking. Results indicated that chlortetracycline quenches BSA fluorescence mainly by a static quenching mechanism. The quenching constants (K_SV) were obtained as 5.64 × 10⁴, 4.49 × 10^4/, and 3.44 × 10^4/ M_?1 at 283, 295, and 307 K, respectively. The thermodynamic parameters of enthalpy change Δ H°, entropy change Δ S°, and free energy change Δ G° were ?5.12 × 10^4/ J mol?1, ?97.6 J mol?1 K?1, and ?2.24 × 10^4/ J mol?1 (295 K), respectively. The association constant (K_A) and the number of binding sites (n) were 9.41 × 10^3/ M?1 and 0.86, respectively. The analysis results suggested that the interaction was spontaneous, and van der Waals force and hydrogen‐bonding interactions played key roles in the reaction process. In addition, CD spectra proved secondary structure alteration of BSA in the presence of chlortetracycline. © 2012 Wiley Periodicals, Inc. J Biochem Mol Toxicol 26:331–336, 2012; View this article online at wileyonlinelibrary.com . DOI 10:1002/jbt.21424     

A molecular dynamics analysis of the GCC-box binding domain in ethylene-responsive element binding factors

Ethylene-responsive element (ERE) binding factors is responsible for a consensus nucleotide sequence AGCCGCC (GCC-box) binding in many important process of plant growing through gene regulation and mediating signal transduction pathways in response to environmental stress. The GCC-box binding domain (GBD) as a novel fold for DNA recognition has been analyzed by means of molecular dynamics. The simulations show that the complex of GBD-DNA trajectories show similar fluctuations in the atomic positions as uncomplexed, particularly at three beta strands involving DNA binding. The calculations of entropy also affirm that GBD flexibility is basically similar for two ligation states. Further, the two complexation states present similar patterns of concerted motions, indicating that the bound DNA cannot alter GBD flexibility. It is inferred that the flexibility of GBD molecule is independent of its ligation state. So in the protein-DNA recognition, the GBD cannot be easily induced while DNA shows better flexibility. Comparison between simulations of unligated GBD and the complexed GBD (in isolation or DNA-bound) reveals intrinsic flexibilities in some certain parts of the molecule play a key role in DNA recognition. In addition, MD simulation identifies that water molecule may mediate interaction between GBD and DNA.