Using molecular dynamics for the refinement of atomistic models of GPCRs by homology modeling

Despite GPCRs sharing a common seven helix bundle, analysis of the diverse crystallographic structures available reveal specific features that might be relevant for ligand design. Despite the number of crystallographic structures of GPCRs steadily increasing, there are still challenges that hamper the availability of new structures. In the absence of a crystallographic structure, homology modeling remains one of the important techniques for constructing 3D models of proteins. In the present study we investigated the use of molecular dynamics simulations for the refinement of GPCRs models constructed by homology modeling. Specifically, we investigated the relevance of template selection, ligand inclusion as well as the length of the simulation on the quality of the GPCRs models constructed. For this purpose we chose the crystallographic structure of the rat muscarinic M3 receptor as reference and constructed diverse atomistic models by homology modeling, using different templates. Specifically, templates used in the present work include the human muscarinic M2; the more distant human histamine H1 and the even more distant bovine rhodopsin as shown in the GPCRs phylogenetic tree. We also investigated the use or not of a ligand in the refinement process. Hence, we conducted the refinement process of the M3 model using the M2 muscarinic as template with tiotropium or NMS docked in the orthosteric site and compared with the results obtained with a model refined without any ligand bound.     

Investigation of ligand selectivity in CYP3A7 by molecular dynamics simulations

Cytochrome P450 (CYP) 3A7 plays a crucial role in the biotransformation of the metabolized endogenous and exogenous steroids. To compare the metabolic capabilities of CYP3A7–ligands complexes, three endogenous ligands were selected, namely dehydroepiandrosterone (DHEA), estrone, and estradiol. In this study, a three-dimensional model of CYP3A7 was constructed by homology modeling using the crystal structure of CYP3A4 as the template and refined by molecular dynamics simulation (MD). The docking method was adopted, combined with MD simulation and the molecular mechanics generalized born surface area method, to probe the ligand selectivity of CYP3A7. These results demonstrate that DHEA has the highest binding affinity, and the results of the binding free energy were in accordance with the experimental conclusion that estrone is better than estradiol. Moreover, several key residues responsible for substrate specificity were identified on the enzyme. Arg372 may be the most important residue due to the low interaction energies and the existence of hydrogen bond with DHEA throughout simulation. In addition, a cluster of Phe residues provides a hydrophobic environment to stabilize ligands. This study provides insights into the structural features of CYP3A7, which could contribute to further understanding of related protein structures and dynamics.     

Insights into the structure and inhibition of Giardia intestinalis arginine deiminase: homology modeling,docking, and molecular dynamics studies

Giardia intestinalis arginine deiminase (GiADI) is an important metabolic enzyme involved in the energy production and defense of this protozoan parasite. The lack of this enzyme in the human host makes GiADI an attractive target for drug design against G. intestinalis. One approach in the design of inhibitors of GiADI could be computer-assisted studies of its crystal structure, such as docking; however, the required crystallographic structure of the enzyme still remains unresolved. Because of its relevance, in this work, we present a three-dimensional structure of GiADI obtained from its amino acid sequence using the homology modeling approximation. Furthermore, we present an approximation of the most stable dimeric structure of GiADI identified through molecular dynamics simulation studies. An in silico analysis of druggability using the structure of GiADI was carried out in order to know if it is a good target for design and optimization of selective inhibitors. Potential GiADI inhibitors were identified by docking of a set of 3196 commercial and 19 in-house benzimidazole derivatives, and molecular dynamics simulation studies were used to evaluate the stability of the ligand–enzyme complexes.     

Structure prediction and molecular dynamics simulations of a G-protein coupled receptor: human CCR2 receptor

CC chemokine receptor type-2 (CCR2) is a member of G-protein coupled receptors superfamily, expressed on the cell surface of monocytes and macrophages. It binds to the monocyte chemoattractant protein-1, a CC chemokine, produced at the sites of inflammation and infection. A homology model of human CCR2 receptor based on the recently available C-X-C chemokine recepor-4 crystal structure has been reported. Ligand information was used as an essential element in the homology modeling process. Six known CCR2 antagonists were docked into the model using simple and induced fit docking procedure. Docked complexes were then subjected to visual inspection to check their suitability to explain the experimental data obtained from site directed mutagenesis and structure-activity relationship studies. The homology model was refined, validated, and assessed for its performance in docking-based virtual screening on a set of CCR2 antagonists and decoys. The docked complexes of CCR2 with the known antagonists, TAK779, a dual CCR2/CCR5 antagonist, and Teijin-comp1, a CCR2 specific antagonist were subjected to molecular dynamics (MD) simulations, which further validated the binding modes of these antagonists. B-factor analysis of 20?ns MD simulations demonstrated that Cys190 is helpful in providing structural rigidity to the extracellular loop (EL2). Residues important for CCR2 antagonism were recognized using free energy decomposition studies. The acidic residue Glu291 from TM7, a conserved residue in chemokine receptors, is favorable for the binding of Teijin-comp1 with CCR2 by ΔG of ?11.4?kcal/mol. Its contribution arises more from the side chains than the backbone atoms. In addition, Tyr193 from EL2 contributes ?0.9?kcal/mol towards the binding of the CCR2 specific antagonist with the receptor. Here, the homology modeling and subsequent molecular modeling studies proved successful in probing the structure of human CCR2 chemokine receptor for the structure-based virtual screening and predicting the binding modes of CCR2 antagonists.     

Computational insights into the binding pattern of mitochondrial calcium uniporter inhibitor through homology modeling,molecular dynamics simulation,binding free energy prediction and density functional theory calculation

Abstract

The mitochondrial calcium uniporter (MCU) is the critical protein of the inner mitochondrial membrane that is the primary mediator for calcium uptake into the mitochondrial matrix. Herein we built the optimal homology model of human MCU which was refined through all-atom molecular dynamics simulation. Then, the binding mode of known inhibitor was predicted through molecular docking method, along with molecular dynamics simulation and binding free energy calculation to verify the docking result and stability of the protein-inhibitor complex. Finally, density functional theory (DFT) calculation enhanced our understanding of the molecular interaction of MCU inhibitor. Our research would provide a deeper insight into the interactions between human MCU and its inhibitor, which boosts to develop novel therapy against MCU related disease.

Communicated by Ramaswamy H. Sarma     

Force-induced unfolding of human telomeric G-quadruplex: A steered molecular dynamics simulation study

We study the unfolding of a parallel G-quadruplex from human telomeric DNA by mechanical stretching using steered molecular dynamics (MD) simulation. We find that the force curves and unfolding processes strongly depend on the pulling sites. With pulling sites located on the sugar-phosphate backbone, the force-extension curve shows a single peak and the unfolding proceeds sequentially. Pulling sites located on the terminal nucleobases lead to a force-extension curve with two peaks and the unfolding is more cooperative. Simulations of the refolding of partially unfolded quadruplexes show very different behavior for the two different pulling modalities. In particular, starting from an unfolded state prepared by nucleobase pulling leads to a long-lived intermediate state whose existence is also corroborated by the free energy profile computed with the Jarzynski equation. Based on this observation, we propose a novel folding pathway for parallel G-quadruplexes with the human telomere sequence.     

A niosomal bilayer of sorbitan monostearate in complex with flavones: a molecular dynamics simulation study

Bilayers prepared from sorbitan fatty acid esters (Span) have been frequently used for delivery of drugs including flavonoids. We applied molecular dynamics simulation to characterize the structure of a sorbitan monostearate (Span 60) bilayer in complex with three representative flavones, a subclass of flavonoids. At a low concentration, unsubstituted flavone, the most hydrophobic member, was able to flip over and cross the bilayer with a large diffusion coefficient. At a high concentration, it was accumulated at the bilayer center resulting in a phase separation. The leaflets of the bilayer were pushed in the opposite directions increasing the membrane thickness. Order parameter of the stearate chain of Span 60 was not affected significantly by unsubstituted flavone. In contrast, chrysin with hydroxylated ring A was lined up with the acyl chains of Span 60 with its hydroxyl group facing the membrane surface. Neither flipping nor transbilayer movement were allowed. Diffusion coefficient was only 15–25% of that of unsubstituted flavone and order parameter decreased with the concentration of chrysin. Luteolin, the most hydroxylated member, interacted mainly with the headgroup of Span 60 and assumed many different orientations without crossing the bilayer. Unlike chrysin and unsubstituted flavone the bilayer integrity was disrupted at 50?mol% luteolin. These behaviors and structures of flavones in a Span 60 bilayer can be accounted for by their hydrophobicity and sites of hydroxylation.     

Homology modeling and molecular dynamics simulation of odonthubuthus doriae (Od1) scorpion toxin in comparison to the BmK M1

All of the α-subgroups share similarity in their sequence and structure but different in the toxicity to various voltage-gated sodium channels (VGSCs). We modeled the first 3D structural model of the Od1 based on BmK M1 using homology modeling. The reliability of model for more investigation and compare to BmK M1 has been examined and confirmed. Then the model structure is further refined by energy minimization and molecular dynamics methods. The purpose of this modeling and simulation is comparison toxicity of two mentioned toxins by investigation structural feature of functional regions including core domain, 5-turn and C-terminal which make NC domain. In the one hand, it is intriguing that Od1 in comparison to BmK M1 shows same solvent accessible surface area (SASA) in 5-turn region but a little more exposed and feasibility (more SASA) in C-terminal region and key functional residues of C-terminal such as positive residues Arg58, lys62 and Arg (His)64. These data suggested that Od1 has similarity with BmK M1 but has more toxicity to sodium channel. In the other hand 5-turn proximity of C-terminal to 5-turn in BmK M1with cis peptide bond is less than Od1 without cis peptide bond which is a confirmation with experimental data about BmK M1.A better understanding of the 3-D structure of Od1and comparison to BmK M1 will be helpful for more investigation of functional characters action of natural toxins with a specialized role for VGSCs.     

A molecular dynamics approach for the generation of complete protein structures from limited coordinate data

Generation of full protein coordinates from limited information, e.g., the Cα coordinates, is an important step in protein homology modeling and structure determination, and molecular dynamics (MD) simulations may prove to be important in this task. We describe a new method, in which the protein backbone is built quickly in a rather crude way and then refined by minimization techniques. Subsequently, the side chains are positioned using extensive MD calculations. The method is tested on two proteins, and results compared to proteins constructed using two other MD-based methods. In the first method, we supplemented an existing backbone building method with a new procedure to add side chains. The second one largely consists of available methodology. The constructed proteins are compared to the corresponding X-ray structures, which became available during this study, and they are in good agreement (backbone RMS values of 0.5–0.7 Å, and all-atom RMS values of 1.5–1.9 Å). This comparative study indicates that extensive MD simulations are able, to some extent, to generate details of the native protein structure, and may contribute to the development of a standardized methodology to predict reliably (parts of) protein structures when only partial coordinate data are available. © 1994 John Wiley & Sons, Inc.     

Atomic insight into designed carbamate-based derivatives as acetylcholine esterase (AChE) inhibitors: a computational study by multiple molecular docking and molecular dynamics simulation

Over 100 variants have been designed and studied, using multiple docking methods such as Autodock Vina, ArgusLab, Molegro Virtual Docker, and Hex-Cuda, to study the effect of alteration in the structure of carbamate-based acetylcholyne esterase (AChE) inhibitors. Sixteen selected systems were then subjected to 14 ns molecular dynamics (MD) simulations. Results from all the docking methods are in agreement. Variants that involved biphenyl substituents possess the most negative binding energies in the ?37.64 to ?39.31 kJ mol^?1 range due to their π–π interactions with AChE aromatic residues. The root mean square deviation values showed that all of these components achieved equilibration after 6 ns. Gyration radius (R_g) and solvent accessibility surface area were calculated to further investigate the AChE conformational changes in the presence of these components. MD simulation results suggested that these components might interact with AChE, possibly with no major changes in AChE secondary and tertiary structures.     

Affinity maturation of a novel antagonistic human monoclonal antibody with a long VH CDR3 targeting the Class A GPCR formyl-peptide receptor 1

Therapeutic monoclonal antibodies targeting G-protein-coupled receptors (GPCRs) are desirable for intervention in a wide range of disease processes. The discovery of such antibodies is challenging due to a lack of stability of many GPCRs as purified proteins. We describe here the generation of Fpro0165, a human anti-formyl peptide receptor 1 (FPR1) antibody generated by variable domain engineering of an antibody derived by immunization of transgenic mice expressing human variable region genes. Antibody isolation and subsequent engineering of affinity, potency and species cross-reactivity using phage display were achieved using FPR1 expressed on HEK cells for immunization and selection, along with calcium release cellular assays for antibody screening. Fpro0165 shows full neutralization of formyl peptide-mediated activation of primary human neutrophils. A crystal structure of the Fpro0165 Fab shows a long, protruding V_H CDR3 of 24 amino acids and in silico docking with a homology model of FPR1 suggests that this long V_H CDR3 is critical to the predicted binding mode of the antibody. Antibody mutation studies identify the apex of the long V_H CDR3 as key to mediating the species cross-reactivity profile of the antibody. This study illustrates an approach for antibody discovery and affinity engineering to typically intractable membrane proteins such as GPCRs.     

Determination of potential selective inhibitors for ROCKI and ROCKII isoforms with molecular modeling techniques: structure based docking,ADMET and molecular dynamics simulation

Rho-associated protein kinases (ROCKs) are a member of the serine/threonine protein kinase family and potential therapeutic target for various diseases. This enzyme has two isoforms, Rho-associated protein kinase I (ROCKI) and Rho-associated protein kinase II (ROCKII). They share an overall 65% homology in all amino acid sequence and 92% homology in kinase domains. Since, the kinase domains of ROCKI and ROCKII are highly conserved and similar, the discovery and design of isoform-selective inhibitors are more challenging. Thus, most currently available agents that is against ROCKs exhibit low selectivity and severe side effects. Therefore, this study aimed to elucidate the interaction of compounds that indicated high potential in experimental studies against ROCKI and ROCKII enzymes in the molecular level with molecular modeling techniques. Firstly, we determined the interaction property of catalytic sites of the ROCKs by analyzing with molecular docking. Based on these results, the best ligands (50 compounds) corresponding to experimental studies were selected, and then absorption, distribution, metabolism and excretion – toxicity (ADMET) analysis of these compounds were implemented. According to these study results, the compound 40 for ROCKI and the compound 50 for ROCKII were identified as selective and highly potent inhibitors. And finally, molecular dynamics (MD) simulations were performed for the stability of ROCKs with identified compounds. In the light of this study, it will be possible to treat diseases that ROCKs have a role by developing more effective and specific ROCK inhibitors.

Communicated by Ramaswamy H. Sarma     

Crystallographic structure versus homology model: a case study of molecular dynamics simulation of human and zebrafish histone deacetylase 10

Abstract

Histone deacetylase (HDAC) 10 has been implicated in the pathology of various cancers and neurodegenerative disorders, making the discovery of novel inhibitors of the isoform an important endeavor. However, the unavailability of crystallographic structure of human HDAC10 (hHDAC10) hinders structure-based drug design effort. Previously, we reported the homology modeled structure of human HDAC10 built using the crystallographic structure of Danio rerio (zebrafish) HDAC10 (zHDAC10) (Protein Data Bank (PDB) ID; 5TD7, released on 24 May 2017) as a template. Here, in continuation with our study, both hHDAC10 and zHDAC10, and their respective complexes with trichostatin A (TSA), quisinostat, and the native ligand (in 5TD7), 7-[(3-aminopropyl)amino]-1,1,1-trifluoroheptane-2,2-diol (PDB ID; FKS) were submitted to 100?ns-long unrestrained molecular dynamics (MD) simulations. Comparative analyses of the MD trajectories revealed that zHDAC10 and its complexes displayed higher stability than hHDAC10 and its corresponding complexes over time. Nonetheless, docking of active and inactive set molecules revealed that more reliable conformations of hHDAC10 could be obtained at an extended time period. This study may shed more light on the reliability of hHDAC10 modeled structure for use in selective inhibitor design.

Communicated by Ramaswamy H. Sarma     

Effect of gold nanoparticle on stability of the DNA molecule: A study of molecular dynamics simulation

An understanding of the mechanism of DNA interactions with gold nanoparticles is useful in today medicine applications. We have performed a molecular dynamics simulation on a B-DNA duplex (CCTCAGGCCTCC) in the vicinity of a gold nanoparticle with a truncated octahedron structure composed of 201 gold atoms (diameter ～1.8 nm) to investigate gold nanoparticle (GNP) effects on the stability of DNA. During simulation, the nanoparticle is closed to DNA and phosphate groups direct the particles into the major grooves of the DNA molecule. Because of peeling and untwisting states that are occur at end of DNA, the nucleotide base lies flat on the surface of GNP. The configuration entropy is estimated using the covariance matrix of atom-positional fluctuations for different bases. The results show that when a gold nanoparticle has interaction with DNA, entropy increases. The results of conformational energy and the hydrogen bond numbers for DNA indicated that DNA becomes unstable in the vicinity of a gold nanoparticle. The radial distribution function was calculated for water hydrogen–phosphate oxygen pairs. Almost for all nucleotide, the presence of a nanoparticle around DNA caused water molecules to be released from the DNA duplex and cations were close to the DNA.     

Effects of cholesterol concentration on the interaction of cytarabine with lipid membranes: a molecular dynamics simulation study

Liposomal cytarabine, DepoCyt, is a chemotherapy agent which is used in cancer treatment. This form of cytarabine has more efficacy and fewer side effects relative to the other forms. Since DepoCyt contains the cytarabine encapsulated within phosphatidylcholine and the sterol molecules, we modeled dioleoylphosphatidylcholine (DOPC)/cholesterol bilayer membrane as a carrier for cytarabine to study drug–bilayer interactions. For this purpose, we performed a series of united-atom molecular dynamics (MD) simulations for 25?ns to investigate the interactions between cytarabine and cholesterol-containing DOPC lipid bilayers. Only the uncharged form of cytarabine molecule was investigated. In this study, different levels of the cholesterol content (0, 20, and 40%) were used. MD simulations allowed us to determine dynamical and structural properties of the bilayer membrane and to estimate the preferred location and orientation of the cytarabine molecule inside the bilayer membrane. Properties such as membrane thickness, area per lipid, diffusion coefficient, mass density, bilayer packing, order parameters, and intermolecular interactions were examined. The results show that by increasing the cholesterol concentration in the lipid bilayers, the bilayer thickness increases and area per lipid decreases. Moreover, in accordance with the experiments, our calculations show that cholesterol molecules have ordering effect on the hydrocarbon acyl chains. Furthermore, the cytarabine molecule preferentially occupies the polar region of the lipid head groups to form specific interactions (hydrogen bonds). Our results fully support the experimental data. Our finding about drug–bilayer interaction is crucial for the liposomal drug design.     

Differential dynamics of the serotonin1A receptor in membrane bilayers of varying cholesterol content revealed by all atom molecular dynamics simulation

Abstract

The serotonin_1A receptor belongs to the superfamily of G protein-coupled receptors (GPCRs) and is a potential drug target in neuropsychiatric disorders. The receptor has been shown to require membrane cholesterol for its organization, dynamics and function. Although recent work suggests a close interaction of cholesterol with the receptor, the structural integrity of the serotonin_1A receptor in the presence of cholesterol has not been explored. In this work, we have carried out all atom molecular dynamics simulations, totaling to 3?μs, to analyze the effect of cholesterol on the structure and dynamics of the serotonin_1A receptor. Our results show that the presence of physiologically relevant concentration of membrane cholesterol alters conformational dynamics of the serotonin_1A receptor and, on an average lowers conformational fluctuations. Our results show that, in general, transmembrane helix VII is most affected by the absence of membrane cholesterol. These results are in overall agreement with experimental data showing enhancement of GPCR stability in the presence of membrane cholesterol. Our results constitute a molecular level understanding of GPCR-cholesterol interaction, and represent an important step in our overall understanding of GPCR function in health and disease.     

Mechanical properties of armchair silicene nanoribbons with edge cracks: a molecular dynamics study

Silicene has been proven to be a promising material with attractive electronic properties. During the synthesis of silicene, structural defects such as edge crack are likely to be generated and such defects in silicene have impacts on its properties. Herein, molecular dynamics simulations were performed to investigate the mechanical properties of the armchair silicene nanoribbons (ASiNRs) with edge cracks. Our results showed that the mechanical properties of the ASiNRs decrease because of the existence of edge crack. Both the pristine ASiNRs and the ASiNRs with edge cracks show brittle fracture behaviours. The crack length plays an important role in determining the critical strain and fracture strength of the ASiNRs. Moreover, we investigated the effects of strain rate and temperature on the mechanical properties of the ASiNRs with edge cracks. We observed that the increasing strain rate increases the critical strain and fracture strength while decreasing the Young’s modulus. Low-strain rates also changes the expanded directions of cracks in the ASiNRs. We also found that the increasing temperature could significantly decrease the mechanical properties of the ASiNRs with edge cracks.     

A molecular dynamics study of liquid aliphatic alcohols: simulation of density and self-diffusion coefficient using a modified OPLS force field

A set of 13 aliphatic alcohols was modelled by molecular dynamics simulations at temperatures from 288 to 338 K using the optimised potential for liquid simulations (OPLS) united-atom force field, the OPLS all-atom force field and the OPLS all-atom force field with modified partial charges of the hydroxyl group. The set includes primary and secondary alcohols, and mono-, di- and trialcohols, and covers a broad range of polarities from log P = ? 0.74 (methanol) to log P = 2.9 (octanol). The density, the radial distribution function, the self-diffusion coefficient and the dielectric constant were evaluated. A long equilibration time of at least 50 ns and a large size of the molecular system of more than 75,000 atoms were used. Except for glycerol, the OPLS all-atom force field reliably reproduced the experimentally determined density with deviations of less than 4% over the whole temperature range. In contrast, the modelled self-diffusion coefficient deviated from its experimental value by up to 55%. To modify the force field, the partial charges of the hydroxyl group were varied by up to 3%. Using the modified OPLS force field, the deviation of the self-diffusion coefficients from their experimental values decreased to less than 19%, while the densities changed by less than 1%.     

Exploring the mechanism of a regulatory SNP of KLK3 by molecular dynamics simulation

The SNP -158G>A of KLK3 has been validated as a regulatory SNP (rSNP) by molecular biology assays, but the mechanism of how it affects the binding of an androgen receptor (AR) homodimer with DNA is unclear. In the current study, molecular dynamics simulation was adopted to explain its inner cause. Based on a recent review), three types of intermolecular forces were analyzed, and the differences among them were compared between complexes containing -158 A:T and -158 G:C. Extra hydrophobic contacts caused by the methyl group on the mutated thymine were the most crucial factor to the regulatory effect of this rSNP. Further analysis concerning the relative motion of the two recognition helixes of the AR homodimer indicated that the hydrophobic interactions between the recognition helix B and the major groove containing -158 A:T changed that helix’s motion greatly from swaying in a plane at free state to vibrating slightly around an equilibrium position. A relatively full explanation on the occurrence of rSNP -158G>A is presented here.     

New insights into the molecular mechanism of methanol-induced inactivation of Thermomyces lanuginosus lipase: a molecular dynamics simulation study

Methanol intolerance of lipase is a major limitation in lipase-catalysed methanolysis reactions. In this study, to understand the molecular mechanism of methanol-induced inactivation of lipases, we performed molecular dynamics (MD) simulations of Thermomyces lanuginosus lipase (TLL) in water and methanol and compared the observed structural and dynamic properties. The solvent accessibility analysis showed that in methanol, polar residues tended to be buried away from the solvent while non-polar residues tended to be more solvent-exposed in comparison to those in water. Moreover, we observed that in methanol, the van der Waals packing of the core residues in two hydrophobic regions of TLL became weak. Additionally, the catalytically relevant hydrogen bond between Asp²⁰¹ OD2 and His²⁵⁸ ND1 in the active site was broken when enzyme was solvated in methanol. This may affect the stability of the tetrahedral intermediates in the catalytic cycle of TLL. Furthermore, compared to in water, some enzyme surface residues displayed enhanced movement in methanol with higher Cα root-mean-square atomic positional fluctuation values. One of such methanol-affecting surface residues (Ile²⁴¹) was chosen for mutation, and MD simulation of the I241E mutant in methanol was conducted. The structural analysis of the mutant showed that replacing a non-polar surface residue with an acidic one at position 241 contributed to the stabilisation of enzyme structure in methanol. Ultimately, these results, while providing molecular-level insights into the destabilising effect of methanol on TLL, highlight the importance of surface residue redesign to improve the stability of lipases in methanol environments.