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

Exhaustively Exploring the Prevalent Interaction Pathways of Ligands Targeting the Ligand-Binding Pocket of Farnesoid X Receptor via Combined Enhanced Sampling

It is well-known that the potency of a drug is heavily associated with its kinetic and thermodynamic properties with the target. Nuclear receptors (NRs), as an important target family, play important roles in regulating a variety of physiological processes in vivo. However, it is hard to understand the drug-NR interaction process because of the closed structure of the ligand-binding domain (LBD) of the NR proteins, which apparently hinders the rational design of drugs with controllable kinetic properties. Therefore, understanding the underlying mechanism of the ligand-NR interaction process seems necessary to help NR drug design. However, it is usually difficult for experimental approaches to interpret the kinetic process of drug-target interactions. Therefore, in silico methods were utilized to explore the optimal binding/dissociation pathways of the NR ligands. Specifically, farnesoid X receptor (FXR) is considered here as the target system since it has been an important target for the treatment of bile acid metabolism-associated diseases, and a series of structures cocrystallized with diverse scaffold ligands were resolved. By using random acceleration molecular dynamics (RAMD) simulation and umbrella sampling (US), 5 main dissociation pathways (pathways I-V) were identified in 11 representative FXR ligands, with most of them (9/11) preferring to go through Pathway III and the remaining two favoring escaping from Pathway I and IV. Furthermore, key residues functioning in the three main dissociation pathways were revealed by the kinetic residue energy analysis (KREA) based on the US trajectories, which may serve as road-marker residues for rapid identification of the (un)binding pathways of FXR ligands. Moreover, the preferred pathways explored by RAMD simulations are in good agreement with the minimum free energy path identified by the US simulations with the Pearson R = 0.76 between the predicted binding affinity and the experimental data, suggesting that RAMD is suitable for applying in large-scale (un)binding-pathway exploration in the case of ligands with obscure binding tunnels to the target.

The prediction of protein-ligand unbinding for modern drug discovery

INTRODUCTION

Drug-target thermodynamic and kinetic information have perennially important roles in drug design. The prediction of protein-ligand unbinding, which can provide important kinetic information, in experiments continues to face great challenges. Uncovering protein-ligand unbinding through molecular dynamics simulations has become efficient and inexpensive with the progress and enhancement of computing power and sampling methods.

AREAS COVERED

In this review, various sampling methods for protein-ligand unbinding and their basic principles are firstly briefly introduced. Then, their applications in predicting aspects of protein-ligand unbinding, including unbinding pathways, dissociation rate constants, residence time and binding affinity, are discussed.

EXPERT OPINION

Although various sampling methods have been successfully applied in numerous systems, they still have shortcomings and deficiencies. Most enhanced sampling methods require researchers to possess a wealth of prior knowledge of collective variables or reaction coordinates. In addition, most systems studied at present are relatively simple, and the study of complex systems in real drug research remains greatly challenging. Through the combination of machine learning and enhanced sampling methods, prediction accuracy can be further improved, and some problems encountered in complex systems also may be solved.

Mechanism of Vitamin D Receptor Ligand-Binding Domain Regulation Studied by gREST Simulations

The vitamin D receptor ligand-binding domain (VDR-LBD) undergoes conformational changes upon ligand binding. In this nuclear receptor family, agonistic or antagonistic activities are controlled by the conformation of the helix (H)12. However, all crystal structures of VDR-LBD reported to date correspond to the active H12 conformation, regardless of agonist/antagonist binding. To understand the mechanism of VDR-LBD regulation structurally, conformational samplings of agonist- and antagonist-bound rat VDR-LBD were performed using the generalized replica exchange with solute tempering (gREST) method. The gREST simulations demonstrated different structural responses of rat VDR-LBD to agonist or antagonist binding, whereas in conventional molecular dynamics simulations, the conformation was the same as that of the crystal structures, regardless of agonist/antagonist binding. In the gREST simulations, a spontaneous conformational change of H12 was observed only for the antagonist complex. The different responses to agonist/antagonist binding were attributed to hydrophobic core formation at the ligand-binding pocket and cooperative rearrangements of H11. The gREST method can be applied to the examination of structure-activity relationships for multiple VDR-LBD ligands.

CaverDock: a molecular docking-based tool to analyse ligand transport through protein tunnels and channels

MOTIVATION

Protein tunnels and channels are key transport pathways that allow ligands to pass between proteins' external and internal environments. These functionally important structural features warrant detailed attention. It is difficult to study the ligand binding and unbinding processes experimentally, while molecular dynamics simulations can be time-consuming and computationally demanding.

RESULTS

CaverDock is a new software tool for analysing the ligand passage through the biomolecules. The method uses the optimized docking algorithm of AutoDock Vina for ligand placement docking and implements a parallel heuristic algorithm to search the space of possible trajectories. The duration of the simulations takes from minutes to a few hours. Here we describe the implementation of the method and demonstrate CaverDock's usability by: (i) comparison of the results with other available tools, (ii) determination of the robustness with large ensembles of ligands and (iii) the analysis and comparison of the ligand trajectories in engineered tunnels. Thorough testing confirms that CaverDock is applicable for the fast analysis of ligand binding and unbinding in fundamental enzymology and protein engineering.

AVAILABILITY AND IMPLEMENTATION

User guide and binaries for Ubuntu are freely available for non-commercial use at https://loschmidt.chemi.muni.cz/caverdock/. The web implementation is available at https://loschmidt.chemi.muni.cz/caverweb/. The source code is available upon request.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Ligand Pathways in Nuclear Receptors

Nuclear receptors (NRs) are ligand-inducible transcription factors that play an essential role in a multitude of physiological processes as well as diseases, rendering them attractive drug targets. Crystal structures revealed the binding site of NRs to be buried in the core of the protein, with no obvious route for ligands to access this cavity. The process of ligand binding is known to be an often-neglected contributor to the efficacy of drug candidates and is thought to influence the selectivity and specificity of NRs. While experimental methods generally fail to highlight the dynamic processes of ligand access or egress on the atomistic scale, computational methods have provided fundamental insight into the pathways connecting the buried binding pocket to the surrounding environment. Methods based on molecular dynamics (MD) and Monte Carlo simulations have been applied to identify pathways and quantify their capability to transport ligands. Here, we systematically review findings of more than 20 years of research in the field, including the applied methodology and controversies. Further, we establish a unified nomenclature to describe the pathways with respect to their location relative to protein secondary structure elements and summarize findings relevant to drug design. Lastly, we discuss the effect of NR interaction partners such as coactivators and corepressors, as well as mutations on the pathways.

Insight into the delivery channel and selectivity of multiple binding sites in bovine serum albumin towards naphthalimide-polyamine derivatives

Naphthalimide derivatives are types of small-molecule anticancer drug candidates; however, their negative factors and potential side effects make their application limited. The pharmacophores select a direct access into the tumor cells as the first choice; this can reduce the side effect of the anti-cancer drugs on the normal cells. Herein, the delivery and binding of the naphthalimide-polyamine complex assisted by the bovine serum albumin (BSA) protein have been studied by combining several molecular dynamic simulations. The plausible transportation channels and the most favorable pathways for the delivery of the naphthalimide-polyamine complex to two drug sites (DSI and DSII), their thermodynamic and dynamic properties and the mechanisms have been discussed in detail. The residues His287 and Phe394 acted as guards in the DSI and DSII, respectively, which played a gating-switch role by flipping the ring from open to close during the compound delivery. The binding mode, binding energy and substituent effects have been also identified. The two drug sites have different preferences towards the compound with the electron-withdrawing and electron-donating substituents, and their strong interactions are more sensitive to the number of the substituent groups. The naphthalimide-polyamine complexes are more likely to choose DSI, both thermodynamically and dynamically, as compared to DSII. This selective specificity of these two drug sites manipulated by the electron-withdrawing and electron-donating substituents is quite promising for the design of new naphthalimide drugs.

Insights into product release dynamics through structural analyses of thymidylate kinase

Several studies on enzyme catalysis have pointed out that the product release event could be a rate limiting step. In this study, we have compared the release event of two products, Adenosine di-phosphate (ADP) and Thymidine di-phosphate (TDP) from the active-site of human and Thermus thermophilus thymidine mono-phosphate kinase (TMPK), referred to as hTMPK and ttTMPK, respectively. TMPK catalyses the conversion of Thymidine mono-phosphate (TMP) to TDP using ATP as phosphoryl donor in the presence of Mg²⁺ ion. Most of the earlier studies on this enzyme have focused on understanding substrate binding and catalysis, but the critical product release event remains elusive. Competitive binding experiments of the substrates and the products using ttTMPK apo crystals have indicated that the substrate (TMP) can replace the bound product (TDP), even in the presence of an ADP molecule. Further, the existing random accelerated molecular dynamics (RAMD) simulation program was modified to study the release of both the products simultaneously from the active site. The RAMD simulations on product-bound structures of both ttTMPK and hTMPK, revealed that while several exit patterns of the products are permissible, the sequential exit mode is the most preferred pattern for both ttTMPK and hTMPK enzymes. Additionally, the product release from the hTMPK was found to be faster and more directional as compared to ttTMPK. Structural investigation revealed that the critical changes in the residue composition in the LID-region of ttTMPK and hTMPK have an effect on the product release and can be attributed to the observed differences during product release event. Understanding of these dissimilarities is of considerable utility in designing potent inhibitors or prodrugs that can distinguish between eukaryotic and prokaryotic homologues of thymidylate kinase.

Estimation of Drug-Target Residence Times by τ-Random Acceleration Molecular Dynamics Simulations

Drug-target residence time (τ), one of the main determinants of drug efficacy, remains highly challenging to predict computationally and, therefore, is usually not considered in the early stages of drug design. Here, we present an efficient computational method, τ-random acceleration molecular dynamics (τRAMD), for the ranking of drug candidates by their residence time and obtaining insights into ligand-target dissociation mechanisms. We assessed τRAMD on a data set of 70 diverse drug-like ligands of the N-terminal domain of HSP90α, a pharmaceutically important target with a highly flexible binding site, obtaining computed relative residence times with an accuracy of about 2.3τ for 78% of the compounds and less than 2.0τ within congeneric series. Analysis of dissociation trajectories reveals features that affect ligand unbinding rates, including transient polar interactions and steric hindrance. These results suggest that τRAMD will be widely applicable as a computationally efficient aid to improving drug residence times during lead optimization.

QM/MM and MM MD Simulations on the Pyrimidine-Specific Nucleoside Hydrolase: A Comprehensive Understanding of Enzymatic Hydrolysis of Uridine

The pyrimidine-specific nucleoside hydrolase Yeik (CU-NH) from Escherichia coli cleaves the N-glycosidic bond of uridine and cytidine with a 10²-10⁴-fold faster rate than that of purine nucleoside substrates, such as inosine. Such a remarkable substrate specificity and the plausible hydrolytic mechanisms of uridine have been explored by using QM/MM and MM MD simulations. The present calculations show that the relatively stronger hydrogen-bond interactions between uridine and the active-site residues Gln227 and Tyr231 in CU-NH play an important role in enhancing the substrate binding and thus promoting the N-glycosidic bond cleavage, in comparison with inosine. The estimated energy barrier of 30 kcal/mol for the hydrolysis of inosine is much higher than 22 kcal/mol for uridine. Extensive MM MD simulations on the transportation of substrates to the active site of CU-NH indicate that the uridine binding is exothermic by ∼23 kcal/mol, more remarkable than inosine (∼12 kcal/mol). All of these arise from the noncovalent interactions between the substrate and the active site featured in CU-NH, which account for the substrate specificity. Quite differing from other nucleoside hydrolases, here the enzymatic N-glycosidic bond cleavage of uridine is less influenced by its protonation.

Estimation of kinetic and thermodynamic ligand-binding parameters using computational strategies

Kinetic and thermodynamic ligand–protein binding parameters are gaining growing importance as key information to consider in drug discovery. The determination of the molecular structures, using particularly x-ray and NMR techniques, is crucial for understanding how a ligand recognizes its target in the final binding complex. However, for a better understanding of the recognition processes, experimental studies of ligand–protein interactions are needed. Even though several techniques can be used to investigate both thermodynamic and kinetic profiles for a ligand–protein complex, these procedures are very often laborious, time consuming and expensive. In the last 10 years, computational approaches have enormous potential in providing insights into each of the above effects and in parsing their contributions to the changes in both kinetic and thermodynamic binding parameters. The main purpose of this review is to summarize the state of the art of computational strategies for estimating the kinetic and thermodynamic parameters of a ligand–protein binding.

Computational study on the unbinding pathways of B-RAF inhibitors and its implication for the difference of residence time: insight from random acceleration and steered molecular dynamics simulations

Random acceleration and steered molecular dynamics simulations reveal the unbinding pathway of B-RAF inhibitors and the difference in the residence time.

A ligand-entry surface of the nuclear receptor superfamily consists of the helix H3 of the ligand-binding domain

We successfully simulated receptor-ligand complex holo-form formation using the human retinoid X receptor-α ligand-binding domain (LBD) and its natural ligand, 9-cis retinoic acid. The success of this simulation was strongly dependent on the findings for an initial structure between the apo-LBD and the ligand as well as the discovery of the driving forces underlying the ligand-trapping and subsequent ligand-induction processes. Here, we would like to propose the "helix H3 three-point initial-binding hypothesis," which was instrumental in simulating the nuclear receptor (NR) superfamily. Using this hypothesis, we also succeeded in simulating holo-form formation of the human retinoic acid receptor-γ LBD and its natural ligand, all-trans retinoic acid. It is hoped that this hypothesis will facilitate novel understanding of both the ligand-trapping mechanism and the simultaneous C-terminal folding process in NR LBDs, as well as provide a new approach to drug design using a structure-based perspective.

Extracting ligands from receptors by reversed targeted molecular dynamics

Short targeted MD trajectories are used to expel ligands from binding sites. The expulsion is governed by a linear increase of the target RMSD value, growing from zero to an arbitrary chosen final RMSD that forces the ligand to a selected distance outside of the receptor. The RMSD lag (i.e., the difference between the imposed and the actual RMSD) can be used to follow barriers encountered by the ligand during its way out of the receptor. The force constant used for the targeted MD can transform the RMSD lag into a strain energy. Integration of the (time-dependent) strain energy over time yields a value with the dimensions of "action" (i.e, energy multiplied by time) and can serve as a measure for the overall effort required to extract the ligand from its binding site. Possibilities to compare (numerically and graphically) the randomly detected exit pathways are discussed. As an example, the method is tested on the exit of bisphenol A from the human estrogen-related receptor [Formula: see text] and of GW0072 from the peroxysome proliferator activated receptor.

Computational insights into inhibitory mechanism of azole compounds against human aromatase

We investigated the inhibitory mechanism of azole aromatase inhibitors. The results showed that letrozole and imazalil prefer different unbinding pathways.

Allosteric mechanisms of nuclear receptors: insights from computational simulations

The traditional structural view of allostery defines this key regulatory mechanism as the ability of one conformational event (allosteric site) to initiate another in a separate location (active site). In recent years computational simulations conducted to understand how this phenomenon occurs in nuclear receptors (NRs) has gained significant traction. These results have yield insights into allosteric changes and communication mechanisms that underpin ligand binding, coactivator binding site formation, post-translational modifications, and oncogenic mutations. Moreover, substantial efforts have been made in understanding the dynamic processes involved in ligand binding and coregulator recruitment to different NR conformations in order to predict cell/tissue-selective pharmacological outcomes of drugs. They also have improved the accuracy of in silico screening protocols so that nowadays they are becoming part of optimisation protocols for novel therapeutics. Here we summarise the important contributions that computational simulations have made towards understanding the structure/function relationships of NRs and how these can be exploited for rational drug design.

Structural insights into the molecular mechanism of vitamin D receptor activation by lithocholic acid involving a new mode of ligand recognition

The vitamin D receptor (VDR), an endocrine nuclear receptor for 1α,25-dihydroxyvitamin D3, acts also as a bile acid sensor by binding lithocholic acid (LCA). The crystal structure of the zebrafish VDR ligand binding domain in complex with LCA and the SRC-2 coactivator peptide reveals the binding of two LCA molecules by VDR. One LCA binds to the canonical ligand-binding pocket, and the second one, which is not fully buried, is anchored to a site located on the VDR surface. Despite the low affinity of the alternative site, the binding of the second molecule promotes stabilization of the active receptor conformation. Biological activity assays, structural analysis, and molecular dynamics simulations indicate that the recognition of two ligand molecules is crucial for VDR agonism by LCA. The unique binding mode of LCA provides clues for the development of new chemical compounds that target alternative binding sites for therapeutic applications.

Dynamics of nuclear receptor Helix-12 switch of transcription activation by modeling time-resolved fluorescence anisotropy decays

Nuclear hormone receptors (NRs) are major targets for pharmaceutical development. Many experiments demonstrate that their C-terminal Helix (H12) is more flexible in the ligand-binding domains (LBDs) without ligand, this increased mobility being correlated with transcription repression and human diseases. Crystal structures have been obtained in which the H12 is extended, suggesting the possibility of large amplitude H12 motions in solution. However, these structures were interpreted as possible crystallographic artifacts, and thus the microscopic nature of H12 movements is not well known. To bridge the gap between experiments and molecular models and provide a definitive picture of H12 motions in solution, extensive molecular dynamics simulations of the peroxisome proliferator-activated receptor-γ LBD, in which the H12 was bound to a fluorescent probe, were performed. A direct comparison of the modeled anisotropy decays to time-resolved fluorescence anisotropy experiments was obtained. It is shown that the decay rates are dependent on the interactions of the probe with the surface of the protein, and display little correlation with the flexibility of the H12. Nevertheless, for the probe to interact with the surface of the LBD, the H12 must be folded over the body of the LBD. Therefore, the molecular mobility of the H12 should preserve the globularity of the LBD, so that ligand binding and dissociation occur by diffusion through the surface of a compact receptor. These results advance the comprehension of both ligand-bound and ligand-free receptor structures in solution, and also guide the interpretation of time-resolved anisotropy decays from a molecular perspective, particularly by the use of simulations.

Identification of a new hormone-binding site on the surface of thyroid hormone receptor

Thyroid hormone receptors (TRs) are members of the nuclear receptor superfamily of ligand-activated transcription factors involved in cell differentiation, growth, and homeostasis. Although X-ray structures of many nuclear receptor ligand-binding domains (LBDs) reveal that the ligand binds within the hydrophobic core of the ligand-binding pocket, a few studies suggest the possibility of ligands binding to other sites. Here, we report a new x-ray crystallographic structure of TR-LBD that shows a second binding site for T3 and T4 located between H9, H10, and H11 of the TRα LBD surface. Statistical multiple sequence analysis, site-directed mutagenesis, and cell transactivation assays indicate that residues of the second binding site could be important for the TR function. We also conducted molecular dynamics simulations to investigate ligand mobility and ligand-protein interaction for T3 and T4 bound to this new TR surface-binding site. Extensive molecular dynamics simulations designed to compute ligand-protein dissociation constant indicate that the binding affinities to this surface site are of the order of the plasma and intracellular concentrations of the thyroid hormones, suggesting that ligands may bind to this new binding site under physiological conditions. Therefore, the second binding site could be useful as a new target site for drug design and could modulate selectively TR functions.

Molecular modeling revealed that ligand dissociation from thyroid hormone receptors is affected by receptor heterodimerization

Numerous ligands bind tightly to thyroid hormone receptors (TRs), and exploring the binding and dissociation of these ligands from TRs will increase our understanding of their mechanisms of action. TRs form transcriptionally active heterodimers with retinoid X receptor (RXR); whether this heterodimerization affects ligand dissociation is poorly understood. To investigate the effects of heterodimerization, classical molecular dynamics (MD) simulations and random acceleration molecular dynamics (RAMD) simulations were performed to probe the dissociation of triiodothyronine (T3) from a TRα-RXR ligand binding domain (LBD) heterodimer and the TRα and TRβ LBDs at the atomic level. Seven (I-VII) dissociation pathways were identified for T3. Heterodimerization inhibited pathway I in the TRα-RXR LBD heterodimer, which may block the proper orientation of the helix 12 (H12), therefore affecting the biological functions of TRs. Upon TR heterodimerization, the second most dominant dissociation pathway switched from pathway IV for TRα LBD to pathway II for TRα-RXR LBD. No significant effects of TR heterodimerization were observed on the dominant dissociation pathway III that was located between H3, the H1-H2 loop and the β-sheet. Our study revealed that TR heterodimerization significantly affects T3 dissociation, which provides important information for the study of other TR ligands.

Agonist and antagonist binding to the nuclear vitamin D receptor: dynamics, mutation effects and functional implications

Purpose

The thermodynamically favored complex between the nuclear vitamin D receptor (VDR) and 1α,25(OH)₂-vitamin D₃ (1,25D3) triggers a shift in equilibrium to favor VDR binding to DNA, heterodimerization with the nuclear retinoid x receptor (RXR) and subsequent regulation of gene transcription. The key amino acids and structural requirements governing VDR binding to nuclear coactivators (NCoA) are well defined. Yet very little is understood about the internal changes in amino acid flexibility underpinning the control of ligand affinity, helix 12 conformation and function. Herein, we use molecular dynamics (MD) to study how the backbone and side-chain flexibility of the VDR differs when a) complexed to 1α,25(OH)₂-vitamin D₃ (1,25D3, agonist) and (23S),25-dehydro-1α(OH)-vitamin D₃-26,23-lactone (MK, antagonist); b) residues that form hydrogen bonds with the C25-OH (H305 and H397) of 1,25D3 are mutated to phenylalanine; c) helix 12 conformation is changed and ligand is removed; and d) x-ray water near the C1- and C3-OH groups of 1,25D3 are present or replaced with explicit solvent.

Methods

We performed molecular dynamic simulations on the apo- and holo-VDRs and used T-Analyst to monitor the changes in the backbone and side-chain flexibility of residues that form regions of the VDR ligand binding pocket (LBP), NCoA surface and control helix 12 conformation.

Results

The VDR-1,25D3 and VDR-MK MD simulations demonstrate that 1,25D3 and MK induce highly similar changes in backbone and side-chain flexibility in residues that form the LBP. MK however did increase the backbone and side-chain flexibility of L404 and R274 respectively. MK also induced expansion of the VDR charge clamp (i.e. NCoA surface) and weakened the intramolecular interaction between H305---V418 (helix 12) and TYR401 (helix 11). In VDR_FF, MK induced a generally more rigid LBP and stronger interaction between F397 and F422 than 1,25D3, and reduced the flexibility of the R274 side-chain. Lastly the VDR MD simulations indicate that R274 can sample multiple conformations in the presence of ligand. When the R274 is extended, the β-OH group of 1,25D3 lies proximal to the backbone carbonyl oxygen of R274 and the side-chain forms H-bonds with hinge domain residues. This differs from the x-ray, kinked geometry, where the side-chain forms an H-bond with the 1α-OH group. Furthermore, 1,25D3, but not MK was observed to stabilize the x-ray geometry of R274 during the > 30 ns MD runs.

Conclusions

The MD methodology applied herein provides an in silico foundation to be expanded upon to better understand the intrinsic flexibility of the VDR and better understand key side-chain and backbone movements involved in the bimolecular interaction between the VDR and its’ ligands.

Electronic supplementary material

The online version of this article (doi:10.1186/2193-9616-1-2) contains supplementary material, which is available to authorized users.

On CO egress from, and re-uptake by, the enzyme MauG, as a mimic of the acquisition of oxidizing agents by the pre-MADH_MauG system. A molecular mechanics approach

In this work, by applying a non-deterministic, randomly-oriented minimal force to the dissociated CO ligand of the MauG-CO system, the molecular-dynamics (MD) behavior of this system could be quickly unraveled. It turned out that CO has no marked directional egress from the high-spin c-heme iron distal pocket. Rather, CO is able to exploit all interstices created during the protein fluctuations. Nonetheless, no steady route toward the surrounding solvent was ever observed: CO jumped first into other binding pockets before being able to escape the protein. In a few cases, on hitting the surrounding H(2)O molecules, CO was observed to reverse direction, re-entering the protein. A contention that conformational inversion of the P107 ring provides a gate to the iron ion is not supported by the present simulations.

Unbinding pathways of GW4064 from human farnesoid X receptor as revealed by molecular dynamics simulations

Farnesoid X receptor (FXR, NR1H4) is a member of a nuclear receptor superfamily, which plays important roles in bile acid homeostasis, lipoprotein and glucose metabolism, and hepatic regeneration. GW4064 is a potent and selective FXR agonist and has become a tool compound to probe the physiological functions of FXR. Until now, the mechanism of GW4064 entering and leaving the FXR pocket is still poorly understood. Here, we report a computational study of GW4064 unbinding pathways from FXR by using several molecular dynamics (MD) simulation techniques. Based on the crystal structure of FXR in complex with GW4064, conventional MD was first used to refine the binding and check the stability of GW4064 in the FXR pocket. Random acceleration MD simulations were then performed to explore the possible unbinding pathways of GW4064 from FXR. Four main pathway clusters were found, among which three subpathways, namely Paths 2A, 2B, and 1B, were observed most frequently. Multiple steered MD simulations were further employed to estimate the maximum rupture force and the sum of the forces and to characterize the intermediate states of the ligand unbinding process. By comparing the average force profiles and structural changes, Paths 2A and 2B were identified to be the most favorable unbinding pathways. The former is located between the H1-H2 loop and the H5-H6 loop, and the latter is located in the cleft formed by the H5-H6 loop, H6, and H7. Moreover, the residues lining the pathways were analyzed for their roles in ligand unbinding. Based on our results, the possible structural modification strategies on GW4064 were also proposed.

Putative binding sites, and pathways to them, for amidine and guanidine current inhibitors on acid-sensing ion channels (ASIC). A theoretical approach with hASIC1a homology model

Central inhibition of the acid-sensing hASIC1a channel, acting upstream of the opiate system, might serve to treat any type of pain, avoiding the unwanted addiction problems of the opioid drugs. To this end, inhibition of hASIC1a channel by PcTx1, a peptide from the Trinidad chevron tarantula, is under development. New inhibitors of the hASIC1a channel are also being sought, in the hope of further modulating the activity, from which antiplasmodial amidine and guanidine phenyl drugs have emerged as promising candidates. However, how such current inhibition takes place remains obscure from the molecular point of view, hindering any further progress in developing drugs. Therefore, the nature of the binding sites, and how they are reached by the amidine-guanidine drugs, was investigated here via automated docking and molecular dynamics with hASIC1a homology models. This study has revealed that this ion channel is rich in binding sites, and that flexible drugs, such as nafamostat, may penetrate it in a snake-like elongated conformation. Then, crawling like a snake through temporary holes in the protein, nafamostat either simply flips, or changes to a high-energy folded conformation to become adapted to the shape of the binding site.

Insights from molecular dynamics: the binding site of cocaine in the dopamine transporter and permeation pathways of substrates in the leucine and dopamine transporters

The dopamine transporter (DAT) facilitates the regulation of synaptic neurotransmitter levels. As a target for therapeutic and illicit psycho-stimulant drugs like antidepressants and cocaine, DAT has been studied intensively. Despite a wealth of mutational and physiological data regarding DAT, the structure remains unsolved and details of the transport mechanism, binding sites and conformational changes remain debated. A bacterial homolog of DAT, the leucine transporter (LeuT(Aa)) has been used as a template and framework for modeling and understanding DAT. Free energy profiles obtained from Multi-Configuration Thermodynamic Integration simulations allowed us to correctly identify the primary and secondary binding pockets of LeuT(Aa). A comparison of free energy profiles for dopamine and cocaine in DAT suggests that the binding site of cocaine is located in a secondary pocket, not the primary substrate site. Two recurring primary pathways for intracellular substrate release from the primary pocket are identified in both transporters using the Random Acceleration Molecular Dynamics method. One pathway appears to follow transmembranes (TMs) 1a and 6b while the other pathway follows along TMs 6b and 8. Interestingly, we observe that a single sodium ion is co-transported with leucine during both simulation types.

Steered molecular dynamics simulations of ligand-receptor interaction in lipocalins

Retinol binding protein (RBP) and an engineered lipocalin, DigA16, have been studied using molecular dynamics simulations. Special emphasis has been placed on explaining the ligand-receptor interaction in RBP-retinol and DigA16-digoxigenin complexes, and steered molecular dynamics simulations of 10-20 ns have been carried out for the ligand expulsion process. Digoxigenin is bound deep inside the cavity of DigA16 and forms several stable hydrogen bonds in addition to the hydrophobic van der Waals interaction with the aromatic side-chains. Four crystalline water molecules inside the ligand-binding cavity remain trapped during the simulations. The strongly hydrophobic receptor site of RBP differs considerably from DigA16, and the main source of ligand attraction comes from the phenyl side-chains. The hydrogen bonds between digoxigenin and DigA16 cause the rupture forces on ligand removal in DigA16 and RBP to differ. The mutated DigA16 residues contribute approximately one-half of the digoxigenin interaction energy with DigA16 and, of these, the energetically most important are residues His35, Arg58, Ser87, Tyr88, and Phe114. Potential "sensor loops" were found for both receptors. These are the outlier loops between residues 114-121 and 63-67 for DigA16 and RBP, respectively, and they are located near the entrance of the ligand-binding cavity. Especially, the residues Glu119 (DigA16) and Leu64 (RBP) are critical for sensing. The ligand binding energies have been estimated based on the linear response approximation of binding affinity by using a previous parametrization for retinoids and RBP.

Protein dynamics and ligand migration interplay as studied by computer simulation

Since proteins are dynamic systems in living organisms, the employment of methodologies contemplating this crucial characteristic results fundamental to allow revealing several aspects of their function. In this work, we present results obtained using classical mechanical atomistic simulation tools applied to understand the connection between protein dynamics and ligand migration. Firstly, we will present a review of the different sampling schemes used in the last years to obtain both ligand migration pathways and the thermodynamic information associated with the process. Secondly, we will focus on representative examples in which the schemes previously presented are employed, concerning the following: i) ligand migration, tunnels, and cavities in myoglobin and neuroglobin; ii) ligand migration in truncated hemoglobin members; iii) NO escape and conformational changes in nitrophorins; iv) ligand selectivity in catalase and hydrogenase; and v) larger ligand migration: the P450 and haloalkane dehalogenase cases. This article is part of a Special Issue entitled: Protein Dynamics: Experimental and Computational Approaches.

Structural modeling of high-affinity thyroid receptor-ligand complexes

Understanding the molecular basis of the binding modes of natural and synthetic ligands to nuclear receptors is fundamental to our comprehension of the activation mechanism of this important class of hormone regulated transcription factors and to the development of new ligands. Thyroid hormone receptors (TR) are particularly important targets for pharmaceuticals development because TRs are associated with the regulation of metabolic rates, body weight, and circulating levels of cholesterol and triglycerides in humans. While several high-affinity ligands are known, structural information is only partially available. In this work we obtain structural models of several TR-ligand complexes with unknown structure by docking high affinity ligands to the receptors' ligand binding domain with subsequent relaxation by molecular dynamics simulations. The binding modes of these ligands are discussed providing novel insights into the development of TR ligands. The experimental binding free energies are reasonably well-reproduced from the proposed models using a simple linear interaction energy free-energy calculation scheme.

Forced unbinding of GPR17 ligands from wild type and R255I mutant receptor models through a computational approach

Background

GPR17 is a hybrid G-protein-coupled receptor (GPCR) activated by two unrelated ligand families, extracellular nucleotides and cysteinyl-leukotrienes (cysteinyl-LTs), and involved in brain damage and repair. Its exploitment as a target for novel neuro-reparative strategies depends on the elucidation of the molecular determinants driving binding of purinergic and leukotrienic ligands. Here, we applied docking and molecular dynamics simulations (MD) to analyse the binding and the forced unbinding of two GPR17 ligands (the endogenous purinergic agonist UDP and the leukotriene receptor antagonist pranlukast from both the wild-type (WT) receptor and a mutant model, where a basic residue hypothesized to be crucial for nucleotide binding had been mutated (R255I) to Ile.

Results

MD suggested that GPR17 nucleotide binding pocket is enclosed between the helical bundle and extracellular loop (EL) 2. The driving interaction involves R255 and the UDP phosphate moiety. To support this hypothesis, steered MD experiments showed that the energy required to unbind UDP is higher for the WT receptor than for R255I. Three potential binding sites for pranlukast where instead found and analysed. In one of its preferential docking conformations, pranlukast tetrazole group is close to R255 and phenyl rings are placed into a subpocket highly conserved among GPCRs. Pulling forces developed to break polar and aromatic interactions of pranlukast were comparable. No differences between the WT receptor and the R255I receptor were found for the unbinding of pranlukast.

Conclusions

These data thus suggest that, in contrast to which has been hypothesized for nucleotides, the lack of the R255 residue doesn't affect the binding of pranlukast a crucial role for R255 in binding of nucleotides to GPR17. Aromatic interactions are instead likely to play a predominant role in the recognition of pranlukast, suggesting that two different binding subsites are present on GPR17.

Computational insights into the mechanism of ligand unbinding and selectivity of estrogen receptors

Estrogen receptors (ER) belong to the nuclear receptor superfamily, and two subtypes, ERalpha and ERbeta, have been identified to date. The differentiated functions and receptor expressions of ERalpha and ERbeta made it attracted to discover subtype-specified ligands with high selectivity. However, these two subtypes are highly homologous and only two residues differ in the ligand binding pocket. Therefore, the mechanism of ligand selectivity has become an important issue in searching selective ligands of ER subtypes. In this study, steered molecular dynamics simulations were carried out to investigate the unbinding pathways of two selective ERbeta ligands from the binding pocket of both ERalpha and ERbeta, which demonstrated that the pathway between the H11 helix and the H7 approximately H8 loop was the most probable for ligand escaping. Then potentials of mean force for ligands unbinding along this pathway were calculated in order to gain insights into the molecular basis for energetics of ligand unbinding and find clues of ligand selectivity. The results indicated that His524/475 in ERalpha/ERbeta acted as a "gatekeeper" during the ligand unbinding. Especially, the H7 approximately H8 loop of ERbeta acted as a polar "transmitter" that controlled the ligand unbinding from the binding site and contributed to the ligand selectivity. Finally, the mechanism of ligand selectivity of ER subtypes was discussed from a kinetic perspective and suggestions for improving the ligand selectivity of ERbeta were also presented. These findings could be helpful for rational design of highly selective ERbeta ligands.

Ligand entry and exit pathways in the beta2-adrenergic receptor

The recently determined crystal structure of the human beta(2)-adrenergic (beta(2)AR) G-protein-coupled receptor provides an excellent structural basis for exploring beta(2)AR-ligand binding and dissociation process. Based on this crystal structure, we simulated ligand exit from the beta(2)AR receptor by applying the random acceleration molecular dynamics (RAMD) simulation method. The simulation results showed that the extracellular opening on the receptor surface was the most frequently observed egress point (referred to as pathway A), and a few other pathways through interhelical clefts were also observed with significantly lower frequencies. In the egress trajectories along pathway A, the D192-K305 salt bridge between the extracellular loop 2 (ECL2) and the apex of the transmembrane helix 7 (TM7) was exclusively broken. The spatial occupancy maps of the ligand computed from the 100 RAMD simulation trajectories indicated that the receptor-ligand interactions that restrained the ligand in the binding pocket were the major resistance encountered by the ligand during exit and no second barrier was notable. We next performed RAMD simulations by using a putative ligand-free conformation of the receptor as input structure. This conformation was obtained in a standard molecular dynamics simulation in the absence of the ligand and it differed from the ligand-bound conformation in a hydrophobic patch bridging ECL2 and TM7 due to the rotation of F193 of ECL2. Results from the RAMD simulations with this putative ligand-free conformation suggest that the cleft formed by the hydrophobic bridge, TM2, TM3, and TM7 on the extracellular surface likely serves as a more specific ligand-entry site and the ECL2-TM7 hydrophobic junction can be partially interrupted upon the entry of ligand that pushes F193 to rotate, resulting in a conformation as observed in the ligand-bound crystal structure. These results may help in the design of beta(2)AR-targeting drugs with improved efficacy, as well as in understanding the receptor subtype selectivity of ligand binding in the beta family of the adrenergic receptors that share almost identical ligand-binding pockets, but show notable amino acid sequence divergence in the putative ligand-entry site, including ECL2 and the extracellular end of TM7.

The vitamin D sterol-vitamin D receptor ensemble model offers unique insights into both genomic and rapid-response signaling

Steroid hormones serve as chemical messengers in a wide number of species and target tissues by transmitting signals that result in both genomic and nongenomic responses. Genomic responses are mediated by the formation of a ligand-receptor complex with its cognate steroid hormone nuclear receptor (NR). Nongenomic responses can be mediated at the plasma membrane by a membrane-localized NR. The focus of this Review is on the structural attributes and molecular mechanisms underlying vitamin D sterol (VDS)-vitamin D receptor (VDR) selective and stereospecific regulation of nongenomic and genomic signaling. The VDS-VDR conformational ensemble model describes how VDSs can selectively initiate or block either nongenomic or genomic biological responses by interacting with two VDR ligand-binding pockets, one kinetically favored by 1alpha,25(OH)(2)D(3) (1,25D) and the other thermodynamically favored. We describe the variables that affect the three major elements of the model: the conformational flexibility of the unliganded (apo) protein, the flexibility of the VDS, and the physicochemical selectivity of the VDR genomic pocket (VDR-GP) and alternative pocket (VDR-AP). We also discuss how these three factors collectively provide a rational explanation for the complexities of VDS regulation of cell biology and highlight the current limitations of the model.

Understanding nuclear receptors using computational methods

Nuclear receptors (NRs) are important targets for therapeutic drugs. NRs regulate transcriptional activities through binding to ligands and interacting with several regulating proteins. Computational methods can provide insights into essential ligand-receptor and protein-protein interactions. These in turn have facilitated the discovery of novel agonists and antagonists with high affinity and specificity as well as have aided in the prediction of toxic side effects of drugs by identifying possible off-target interactions. Here, we review the application of computational methods toward several clinically important NRs (with special emphasis on PXR) and discuss their use for screening and predicting the toxic side effects of xenobiotics.