Unbinding of retinoic acid from the retinoic acid receptor by random expulsion molecular dynamics

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.

Membrane-Dependent Binding and Entry Mechanism of Dopamine into Its Receptor

Synaptic neurotransmission has recently been proposed to function via either a membrane-independent or a membrane-dependent mechanism, depending on the neurotransmitter type. In the membrane-dependent mechanism, amphipathic neurotransmitters first partition to the lipid headgroup region and then diffuse along the membrane plane to their membrane-buried receptors. However, to date, this mechanism has not been demonstrated for any neurotransmitter-receptor complex. Here, we combined isothermal calorimetry measurements with a diverse set of molecular dynamics simulation methods to investigate the partitioning of an amphipathic neurotransmitter (dopamine) and the mechanism of its entry into the ligand-binding site. Our results show that the binding of dopamine to its receptor is consistent with the membrane-dependent binding and entry mechanism. Both experimental and simulation results showed that dopamine favors binding to lipid membranes especially in the headgroup region. Moreover, our simulations revealed a ligand-entry pathway from the membrane to the binding site. This pathway passes through a lateral gate between transmembrane alpha-helices 5 and 6 on the membrane-facing side of the protein. All in all, our results demonstrate that dopamine binds to its receptor by a membrane-dependent mechanism, and this is complemented by the more traditional binding mechanism directly through the aqueous phase. The results suggest that the membrane-dependent mechanism is common in other synaptic receptors, too.

Unbinding of Translocator Protein 18 kDa (TSPO) Ligands: From in Vitro Residence Time to in Vivo Efficacy via in Silico Simulations

Translocator protein 18 kDa (TSPO) is a validated pharmacological target for the development of new treatments for neurological disorders. N,N-Dialkyl-2-phenylindol-3-ylglyoxylamides (PIGAs) are effective TSPO modulators and potentially useful therapeutics for the treatment of anxiety, central nervous system pathologies featuring astrocyte loss, and inflammatory-based neuropathologies. For this class of compounds, no correlation exists between the TSPO binding affinity and the corresponding functional efficacy. Rather, their biological effectiveness correlates with the kinetics of the unbinding events and more specifically with the residence time (RT). So far, the structural reasons for the different recorded RT of congeneric PIGAs remain elusive. Here, to understand the different kinetics of PIGAs, their unbinding paths were studied by employing enhanced-sampling molecular dynamics simulations. Results of these studies revealed how subtle structural differences between PIGAs have a substantial effect on the unbinding energetics. In particular, during the egress from the TSPO binding site, slow-dissociating PIGAs find tight interactions with the protein LP1 region thereby determining a long RT. Further support to these findings was achieved by in vivo studies, which demonstrated how the anxiolytic effect observed for the inspected PIGAs correlated with their RT to TSPO.

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.

Computational modeling approaches to quantitative structure-binding kinetics relationships in drug discovery

Simple comparative correlation analyses and quantitative structure-kinetics relationship (QSKR) models highlight the interplay of kinetic rates and binding affinity as an essential feature in drug design and discovery. The choice of the molecular series, and their structural variations, used in QSKR modeling is fundamental to understanding the mechanistic implications of ligand and/or drug-target binding and/or unbinding processes. Here, we discuss the implications of linear correlations between kinetic rates and binding affinity constants and the relevance of the computational approaches to QSKR modeling.

Single or Multiple Access Channels to the CYP450s Active Site? An Answer from Free Energy Simulations of the Human Aromatase Enzyme

Cytochromes P450 (CYP450s), in particular, CYP19A1 and CYP17A1, are key clinical targets of breast and prostate anticancer therapies, critical players in drug metabolism, and their overexpression in tumors is associated with drug resistance. In these enzymes, ligand (substrates, drugs) metabolism occurs in deeply buried active sites accessible only via several grueling channels, whose exact biological role remains unclear. Gaining direct insights on the mechanism by which ligands travel in and out is becoming increasingly important given that channels are involved in the modulation of binding/dissociation kinetics and the specificity of ligands toward a CYP450. This has profound implications for enzymatic efficiency and drug efficacy/toxicity. Here, by applying free energy methods, for a cumulative simulation time of 20 μs, we provide detailed atomistic characterization and free energy profiles of the entry/exit routes preferentially followed by a substrate (androstenedione) and a last-generation inhibitor (letrozole) to/from the catalytic site of CYP19A1 (the human aromatase (HA) enzyme), a key clinical target against breast cancer, studied here as prototypical CYP450. Despite the remarkably different size/shape/hydrophobicity of the ligands, two channels appear accessible to their entrance, while only one exit route appears to be preferential. Our study shows that the preferential paths may be conserved among different CYP450s. Moreover, our results highlight that, at least in the case of HA, ligand channeling is associated with large enzyme structural rearrangements. A wise choice of the computational method and very long simulations are, thus, required to obtain fully converged quantitative free energy profiles, which might be used to design novel biocatalysts or next-generation cytochrome inhibitors with an in silico tuned K_m.

Memetic algorithms for ligand expulsion from protein cavities

Ligand diffusion through a protein interior is a fundamental process governing biological signaling and enzymatic catalysis. A complex topology of channels in proteins leads often to difficulties in modeling ligand escape pathways by classical molecular dynamics simulations. In this paper, two novel memetic methods for searching the exit paths and cavity space exploration are proposed: Memory Enhanced Random Acceleration (MERA) Molecular Dynamics (MD) and Immune Algorithm (IA). In MERA, a pheromone concept is introduced to optimize an expulsion force. In IA, hybrid learning protocols are exploited to predict ligand exit paths. They are tested on three protein channels with increasing complexity: M2 muscarinic G-protein-coupled receptor, enzyme nitrile hydratase, and heme-protein cytochrome P450cam. In these cases, the memetic methods outperform simulated annealing and random acceleration molecular dynamics. The proposed algorithms are general and appropriate in all problems where an accelerated transport of an object through a network of channels is studied.

Exploring the different ligand escape pathways in acylaminoacyl peptidase by random acceleration and steered molecular dynamics simulations

Acylaminoacyl peptidase (APH, EC 3.4.19.1) is a novel class of serine-type protease belonging to the prolyl oligopeptidase (POP) family.

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.

Ligand tunnels in T. brucei and human CYP51: Insights for parasite-specific drug design

BACKGROUND

Cytochrome P450 sterol 14α-demethylase (CYP51) is an essential enzyme for sterol biosynthesis and a target for anti-parasitic drug design. However, the design of parasite-specific drugs that inhibit parasitic CYP51 without severe side effects remains challenging. The active site of CYP51 is situated in the interior of the protein. Here, we characterize the potential ligand egress routes and mechanisms in Trypanosoma brucei and human CYP51 enzymes.

METHODS

We performed Random Acceleration Molecular Dynamics simulations of the egress of four different ligands from the active site of models of soluble and membrane-bound T. brucei CYP51 and of soluble human CYP51.

RESULTS

In the simulations, tunnel 2f, which leads to the membrane, was found to be the predominant ligand egress tunnel for all the ligands studied. Tunnels S, 1 and W, which lead to the cytosol, were also used in T. brucei CYP51, whereas tunnel 1 was the only other tunnel used significantly in human CYP51. The common tunnels found previously in other CYPs were barely used. The ligand egress times were shorter for human than T. brucei CYP51, suggesting lower barriers to ligand passage. Two gating residues, F105 and M460, in T. brucei CYP51 that modulate the opening of tunnels 2f and S were identified.

CONCLUSIONS

Although the main egress tunnel was the same, differences in the tunnel-lining residues, ligand passage and tunnel usage were found between T. brucei and human CYP51s.

GENERAL SIGNIFICANCE

The results provide a basis for the design of selective anti-parasitic agents targeting the ligand tunnels.

Substrate tunnels in enzymes: structure-function relationships and computational methodology

In enzymes, the active site is the location where incoming substrates are chemically converted to products. In some enzymes, this site is deeply buried within the core of the protein, and, in order to access the active site, substrates must pass through the body of the protein via a tunnel. In many systems, these tunnels act as filters and have been found to influence both substrate specificity and catalytic mechanism. Identifying and understanding how these tunnels exert such control has been of growing interest over the past several years because of implications in fields such as protein engineering and drug design. This growing interest has spurred the development of several computational methods to identify and analyze tunnels and how ligands migrate through these tunnels. The goal of this review is to outline how tunnels influence substrate specificity and catalytic efficiency in enzymes with buried active sites and to provide a brief summary of the computational tools used to identify and evaluate these tunnels.

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.

Residue correlation networks in nuclear receptors reflect functional specialization and the formation of the nematode-specific P-box

Background

Nuclear receptors (NRs) are transcription factors which bind small hormones, whose evolutionary history and the presence of different functional surfaces makes them an interesting target for a correlation based analysis.

Results

Correlation analysis of ligand binding domains shows that correlated residue subsets arise from the differences between functional sites in different nuclear receptor subfamilies. For the DNA binding domain, particularly, the analysis shows that the main source of correlation comes from residues that regulate hormone response element specificity, and one of the conserved residue sub-sets arises due to the presence of an unusual sequence for the DNA binding motif known as P-box in nematodes, suggesting the existence of different DBD-DNA specificities in nuclear receptors.

Conclusions

We conclude that DNA specificity and functional surface specialization has independently driven nuclear receptor evolution, and suggest possible binding modes for the class of divergent nuclear receptors in nematodes.

Ligand release mechanisms and channels in histone deacetylases

Exploring the molecular channels of class I histone deacetylases (HDACs) with buried active sites are important to understand their structures and functionalities. In this work, we perform hybrid classical molecular dynamics and random acceleration molecular dynamics simulations to explore the B3N [i.e., (4-(dimethylamino)N-[7(hydroxyamino)-7-oxoheptyle] benzamide)] exit channels in the x-ray crystal structures of HDAC3 and HDAC8 enzymes. Our simulations identify B3N release through four different channels in HDAC3 (denoted as A1, A2, B1, and B2) and HDAC8 (referred as A1, B1, B2, and B3) enzymes, among which egression through channel A1 is more predominant in both the enzymes. This mechanism is similar to ligand release in HDAC1 and HDAC2 described in our previous study and can be the fingerprint ligand release mechanisms in class I HDACs. Ligand release events through B channels, on the other hand, are different among HDAC3 and HDAC8, highlighting the significances of substituted residues in controlling the access to these channels This study reveals a novel aromatic gating mechanism elicited by TYR154-TRP141-TYR111 that controls the B3N access to all the B channels in HDAC8. The TRP141 in HDAC8 is substituted by LEU133 in HDAC3, which do not hinder the access to B channels in HDAC3. However, two hydrogen bonded barricades formed as ARG28-GLY297-GLY295-GLY131 and TRP129-ARG28-ALA130-LEU29-TRP129 obstruct the B3N from exploring the B channels in HDAC3. The structural and dynamical characterizations of molecular channels and ligand unbinding mechanisms reported in this study provide novel structural insights and atomic level perspectives on HDAC3 and HDAC8 enzymes, thereby potentially aiding in the design of more specific HDAC inhibitors.

Modelling and enhanced molecular dynamics to steer structure-based drug discovery

The ever-increasing gap between the availabilities of the genome sequences and the crystal structures of proteins remains one of the significant challenges to the modern drug discovery efforts. The knowledge of structure-dynamics-functionalities of proteins is important in order to understand several key aspects of structure-based drug discovery, such as drug-protein interactions, drug binding and unbinding mechanisms and protein-protein interactions. This review presents a brief overview on the different state of the art computational approaches that are applied for protein structure modelling and molecular dynamics simulations of biological systems. We give an essence of how different enhanced sampling molecular dynamics approaches, together with regular molecular dynamics methods, assist in steering the structure based drug discovery processes.

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.

Molecular mechanism of peroxisome proliferator-activated receptor α activation by WY14643: a new mode of ligand recognition and receptor stabilization

Peroxisome proliferator-activated receptors (PPARs) are members of a superfamily of nuclear transcription factors. They are involved in mediating numerous physiological effects in humans, including glucose and lipid metabolism. PPARα ligands effectively treat dyslipidemia and have significant antiinflammatory and anti-atherosclerotic activities. These effects and their ligand-dependent activity make nuclear receptors obvious targets for drug design. Here, we present the structure of the human PPARα in complex with WY14643, a member of fibrate class of drug, and a widely used PPAR activator. The crystal structure of this complex suggests that WY14643 induces activation of PPARα in an unusual bipartite mechanism involving conventional direct helix 12 stabilization and an alternative mode that involves a second ligand in the pocket. We present structural observations, molecular dynamics and activity assays that support the importance of the second site in WY14643 action. The unique binding mode of WY14643 reveals a new pattern of nuclear receptor ligand recognition and suggests a novel basis for ligand design, offering clues for improving the binding affinity and selectivity of ligand. We show that binding of WY14643 to PPARα was associated with antiinflammatory disease in a human corneal cell model, suggesting possible applications for PPARα ligands.

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.

Alternative pathways for association and dissociation of the calmodulin-binding domain of plasma membrane Ca(2+)-ATPase isoform 4b (PMCA4b)

The calmodulin (CaM)-binding domain of isoform 4b of the plasma membrane Ca(2+) -ATPase (PMCA) pump is represented by peptide C28. CaM binds to either PMCA or C28 by a mechanism in which the primary anchor residue Trp-1093 binds to the C-terminal lobe of the extended CaM molecule, followed by collapse of CaM with the N-terminal lobe binding to the secondary anchor Phe-1110 (Juranic, N., Atanasova, E., Filoteo, A. G., Macura, S., Prendergast, F. G., Penniston, J. T., and Strehler, E. E. (2010) J. Biol. Chem. 285, 4015-4024). This is a relatively rapid reaction, with an apparent half-time of ~1 s. The dissociation of CaM from PMCA4b or C28 is much slower, with an overall half-time of ~10 min. Using targeted molecular dynamics, we now show that dissociation of Ca(2+)-CaM from C28 may occur by a pathway in which Trp-1093, although deeply embedded in a pocket in the C-terminal lobe of CaM, leaves first. The dissociation begins by relatively rapid release of Trp-1093, followed by very slow release of Phe-1110, removal of C28, and return of CaM to its conformation in the free state. Fluorescence measurements and molecular dynamics calculations concur in showing that this alternative path of release of the PMCA4b CaM-binding domain is quite different from that of binding. The intermediate of dissociation with exposed Trp-1093 has a long lifetime (minutes) and may keep the PMCA primed for activation.

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.

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.

An interpretation of positional displacement of the helix12 in nuclear receptors: preexistent swing-up motion triggered by ligand binding

Positional displacement of helix12 (H12) in the estrogen receptor alpha, which belongs to the nuclear receptor (NR) superfamily, is studied by the molecular dynamics (MD) simulation and the linear response theory. Tendency of the H12 to swing up upon ligand binding, which is consistent with X-ray structures and earlier MD simulations, is reproduced by the calculation of the conformational fluctuation in apo state and the response to the external perturbation. Our study thus provides an interpretation of the positional change of the H12 such that it is derived by the preexistent swing-up motion where the ligand binding works only as a trigger. Our finding, which illustrates underlying mechanism of the H12 motion, would contribute to finding a way to regulate the transcriptional activity by synthesized ligands because the transcriptional activity of the NR is governed by the position of the H12.

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.

Ligand unbinding from the estrogen receptor: a computational study of pathways and ligand specificity

The estrogen receptor (ER) belongs to the nuclear receptor superfamily, whose members regulate important cellular events like development and metabolism. The ER functions as a transcription regulator and can be activated on ligand binding. Consequently, ligand binding and unbinding constitute fundamental processes in the regulation of genes. Even though both biochemical and structural data of ER are available, the actual mechanism of the ligand binding/unbinding remains elusive. We have performed computational studies on the unbinding mechanism of ERalpha and ERbeta, in the presence of cofactors and with ligands ranging from agonist to a full antagonist. Our results show that agonists or selective ER modulators can dissociate from the receptor through multiple pathways with minor effect on the receptor structure, whereas an antagonist requires larger conformational changes. Furthermore, a specific receptor/ligand combination can exhibit a pathway preference depending on character and conformation of both parts. Hence, it is possible that the binding/unbinding mechanism can explain ligand subtype specificity and thus have an impact in drug discovery.

Steered molecular dynamics simulations reveal the likelier dissociation pathway of imatinib from its targeting kinases c-Kit and Abl

Development of small molecular kinase inhibitors has recently been the central focus in drug discovery. And type II kinase inhibitors that target inactive conformation of kinases have attracted particular attention since their potency and selectivity are thought to be easier to achieve compared with their counterpart type I inhibitors that target active conformation of kinases. Although mechanisms underlying the interactions between type II inhibitors and their targeting kinases have been widely studied, there are still some challenging problems, for example, how type II inhibitors associate with or dissociate from their targeting kinases. In this investigation, steered molecular dynamics simulations have been carried out to explore the possible dissociation pathways of typical type II inhibitor imatinib from its targeting protein kinases c-Kit and Abl. The simulation results indicate that the most favorable pathway for imatinib dissociation corresponds to the ATP-channel rather than the relatively wider allosteric-pocket-channel, which is mainly due to the different van der Waals interaction that the ligand suffers during dissociation. Nevertheless, the direct reason comes from the fact that the residues composing the ATP-channel are more flexible than that forming the allosteric-pocket-channel. The present investigation suggests that a bulky hydrophobic head is unfavorable, but a large polar tail is allowed for a potent type II inhibitor. The information obtained here can be used to direct the discovery of type II kinase inhibitors.

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.

Efficient table lookup without inverse square roots for calculation of pair wise atomic interactions in classical simulations

A major bottleneck in classical atomistic simulations of biomolecular systems is the calculation of the pair wise nonbonded (Coulomb, van der Waals) interactions. This remains an issue even when methods are used (e.g., lattice summation or spherical cutoffs) in which the number of interactions is reduced from O(N(2)) to O(NlogN) or O(N). The interaction forces and energies can either be calculated directly each time they are needed or retrieved using precomputed values in a lookup table; the choice between direct calculation and table lookup methods depends on the characteristics of the system studied (total number of particles and the number of particle kinds) as well as the hardware used (CPU speed, size and speed of cache, and main memory). A recently developed lookup table code, implemented in portable and easily maintained FORTRAN 95 in the CHARMM program (www.charmm.org), achieves a 1.5- to 2-fold speedup compared with standard calculations using highly optimized FORTRAN code in real molecular dynamics simulations for a wide range of molecular system sizes. No approximations other than the finite resolution of the tables are introduced, and linear interpolation in a table with the relatively modest density of 100 points/A(2) yields the same accuracy as the standard double precision calculations. For proteins in explicit water a less dense table (10 points/A(2)) is 10-20% faster than using the larger table, and only slightly less accurate. The lookup table is even faster than hand coded assembler routines in most cases, mainly due to a significantly smaller operation count inside the inner loop.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Pathways and mechanisms for product release in the engineered haloalkane dehalogenases explored using classical and random acceleration molecular dynamics simulations

Eight mutants of the DhaA haloalkane dehalogenase carrying mutations at the residues lining two tunnels, previously observed by protein X-ray crystallography, were constructed and biochemically characterized. The mutants showed distinct catalytic efficiencies with the halogenated substrate 1,2,3-trichloropropane. Release pathways for the two dehalogenation products, 2,3-dichloropropane-1-ol and the chloride ion, and exchange pathways for water molecules, were studied using classical and random acceleration molecular dynamics simulations. Five different pathways, denoted p1, p2a, p2b, p2c, and p3, were identified. The individual pathways showed differing selectivity for the products: the chloride ion releases solely through p1, whereas the alcohol releases through all five pathways. Water molecules play a crucial role for release of both products by breakage of their hydrogen-bonding interactions with the active-site residues and shielding the charged chloride ion during its passage through a hydrophobic tunnel. Exchange of the chloride ions, the alcohol product, and the waters between the buried active site and the bulk solvent can be realized by three different mechanisms: (i) passage through a permanent tunnel, (ii) passage through a transient tunnel, and (iii) migration through a protein matrix. We demonstrate that the accessibility of the pathways and the mechanisms of ligand exchange were modified by mutations. Insertion of bulky aromatic residues in the tunnel corresponding to pathway p1 leads to reduced accessibility to the ligands and a change in mechanism of opening from permanent to transient. We propose that engineering the accessibility of tunnels and the mechanisms of ligand exchange is a powerful strategy for modification of the functional properties of enzymes with buried active sites.

Molecular dynamics simulation of ligand dissociation from liver fatty acid binding protein

The mechanisms of how ligands enter and leave the binding cavity of fatty acid binding proteins (FABPs) have been a puzzling question over decades. Liver fatty acid binding protein (LFABP) is a unique family member which accommodates two molecules of fatty acids in its cavity and exhibits the capability of interacting with a variety of ligands with different chemical structures and properties. Investigating the ligand dissociation processes of LFABP is thus a quite interesting topic, which however is rather difficult for both experimental approaches and ordinary simulation strategies. In the current study, random expulsion molecular dynamics simulation, which accelerates ligand motions for rapid dissociation, was used to explore the potential egress routes of ligands from LFABP. The results showed that the previously hypothesized “portal region” could be readily used for the dissociation of ligands at both the low affinity site and the high affinity site. Besides, one alternative portal was shown to be highly favorable for ligand egress from the high affinity site and be related to the unique structural feature of LFABP. This result lends strong support to the hypothesis from the previous NMR exchange studies, which in turn indicates an important role for this alternative portal. Another less favored potential portal located near the N-terminal end was also identified. Identification of the dissociation pathways will allow further mechanistic understanding of fatty acid uptake and release by computational and/or experimental techniques.

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.