Multiscale computational study of ligand binding pathways: Case of p38 MAP kinase and its inhibitors

Protein kinases are one of the most important drug targets in the past 10 years. Understanding the inhibitor association processes will profoundly impact new binder designs with preferred binding kinetics. However, after more than a decade of effort, a complete atomistic-level study of kinase inhibitor binding pathways is still lacking. As all kinases share a similar scaffold, we used p38 kinase as a model system to investigate the conformational dynamics and free energy transition of inhibitor binding toward kinases. Two major kinase conformations, Asp-Phe-Gly (DFG)-in and DFG-out, and three types of inhibitors, type I, II, and III, were thoroughly investigated in this work. We performed Brownian dynamics simulations and up to 340 μs Gaussian-accelerated molecular dynamics simulations to capture the inhibitor binding paths and a series of conformational transitions of the p38 kinase from its apo to inhibitor-bound form. Eighteen successful binding trajectories, including all types of inhibitors, are reported herein. Our simulations suggest a mechanism of inhibitor recruitment, a faster ligand association step to a pre-existing DFG-in/DFG-out p38 protein, followed by a slower molecular rearrangement step to adjust the protein-ligand conformation followed by a shift in the energy landscape to reach the final bound state. The ligand association processes also reflect the energetic favor of type I and type II/III inhibitor binding through ATP and allosteric channels, respectively. These different binding routes are directly responsible for the fast (type I binders) and slow (type II/III binders) kinetics of different types of p38 inhibitors. Our findings also echo the recent study of p38 inhibitor dissociation, implying that ligand unbinding could undergo a reverse path of binding, and both processes share similar metastates. This study deepens the understanding of molecular and energetic features of kinase inhibitor-binding processes and will inspire future drug development from a kinetic point of view.     

Switches of hydrogen bonds during ligand–protein association processes determine binding kinetics

Revealing the processes of ligand–protein associations deepens our understanding of molecular recognition and binding kinetics. Hydrogen bonds (H‐bonds) play a crucial role in optimizing ligand–protein interactions and ligand specificity. In addition to the formation of stable H‐bonds in the final bound state, the formation of transient H‐bonds during binding processes contributes binding kinetics that define a ligand as a fast or slow binder, which also affects drug action. However, the effect of forming the transient H‐bonds on the kinetic properties is little understood. Guided by results from coarse‐grained Brownian dynamics simulations, we used classical molecular dynamics simulations in an implicit solvent model and accelerated molecular dynamics simulations in explicit waters to show that the position and distribution of the H‐bond donor or acceptor of a drug result in switching intermolecular and intramolecular H‐bond pairs during ligand recognition processes. We studied two major types of HIV‐1 protease ligands: a fast binder, xk263, and a slow binder, ritonavir. The slow association rate in ritonavir can be attributed to increased flexibility of ritonavir, which yields multistep transitions and stepwise entering patterns and the formation and breaking of complex H‐bond pairs during the binding process. This model suggests the importance of conversions of spatiotemporal H‐bonds during the association of ligands and proteins, which helps in designing inhibitors with preferred binding kinetics. Copyright © 2014 John Wiley & Sons, Ltd.     

SPT and Imaging FCS Provide Complementary Information on the Dynamics of Plasma Membrane Molecules

The dynamics of biomolecules in the plasma membrane is of fundamental importance to understanding cellular processes. Cellular signaling often starts with extracellular ligand binding to a membrane receptor, which then transduces an intracellular signal. Ligand binding and receptor-complex activation often involve a complex rearrangement of proteins in the membrane, which results in changes in diffusion properties. Two widely used methods to characterize biomolecular diffusion are single-particle tracking (SPT) and imaging total internal reflection fluorescence correlation spectroscopy (ITIR-FCS). Here, we compare the results of recovered diffusion coefficients and mean-square displacements of the two methods by simulations of free, domain-confined, or meshwork diffusion. We introduce, to our knowledge, a new method for the determination of confinement radii from ITIR-FCS data. We further establish and demonstrate simultaneous SPT/ITIR-FCS for direct comparison within living cells. Finally, we compare the results obtained by SPT and ITIR-FCS for the receptor tyrosine kinase MET. Our results show that SPT and ITIR-FCS yield complementary information on diffusion properties of biomolecules in cell membranes.     

Towards understanding the mechanisms of molecular recognition by computer simulations of ligand-protein interactions

The thermodynamic and kinetic aspects of molecular recognition for the methotrexate (MTX)-dihydrofolate reductase (DHFR) ligand-protein system are investigated by the binding energy landscape approach. The impact of 'hot' and 'cold' errors in ligand mutations on the thermodynamic stability of the native MTX-DHFR complex is analyzed, and relationships between the molecular recognition mechanism and the degree of ligand optimization are discussed. The nature and relative stability of intermediates and thermodynamic phases on the ligand-protein association pathway are studied, providing new insights into connections between protein folding and molecular recognition mechanisms, and cooperativity of ligand-protein binding. The results of kinetic docking simulations are rationalized based on the thermodynamic properties determined from equilibrium simulations and the shape of the underlying binding energy landscape. We show how evolutionary ligand selection for a receptor active site can produce well-optimized ligand-protein systems such as MTX-DHFR complex with the thermodynamically stable native structure and a direct transition mechanism of binding from unbound conformations to the unique native structure.     

Long-Timescale Molecular-Dynamics Simulations of the Major Urinary Protein Provide Atomistic Interpretations of the Unusual Thermodynamics of Ligand Binding

The mouse major urinary protein (MUP) has proved to be an intriguing test bed for detailed studies on protein-ligand recognition. NMR, calorimetric, and modeling investigations have revealed that the thermodynamics of ligand binding involve a complex interplay between competing enthalpic and entropic terms. We performed six independent, 1.2 μs molecular-dynamics simulations on MUP—three replicates on the apo-protein, and three on the complex with the pheromone isobutylmethoxypyrazine. Our findings provide the most comprehensive picture to date of the structure and dynamics of MUP, and how they are modulated by ligand binding. The mechanical pathways by which amino acid side chains can transmit information regarding ligand binding to surface loops and either increase or decrease their flexibility (entropy-entropy compensation) are identified. Dewetting of the highly hydrophobic binding cavity is confirmed, and the results reveal an aspect of ligand binding that was not observed in earlier, shorter simulations: bound ligand retains extensive rotational freedom. Both of these features have significant implications for interpretations of the entropic component of binding. More generally, these simulations test the ability of current molecular simulation methods to produce a reliable and reproducible picture of protein dynamics on the microsecond timescale.     

Engineering Tocopherol Selectivity in α-TTP: A Combined In Vitro/In Silico Study

We present a combined in vitro/in silico study to determine the molecular origin of the selectivity of -tocopherol transfer protein (-TTP) towards -tocopherol. Molecular dynamics simulations combined to free energy perturbation calculations predict a binding free energy for -tocopherol to -TTP 8.262.13 kcal mol lower than that of -tocopherol. Our calculations show that -tocopherol binds to -TTP in a significantly distorted geometry as compared to that of the natural ligand. Variations in the hydration of the binding pocket and in the protein structure are found as well. We propose a mutation, A156L, which significantly modifies the selectivity properties of -TTP towards the two tocopherols. In particular, our simulations predict that A156L binds preferentially to -tocopherol, with striking structural similarities to the wild-type--tocopherol complex. The affinity properties are confirmed by differential scanning fluorimetry as well as in vitro competitive binding assays. Our data indicate that residue A156 is at a critical position for determination of the selectivity of -TTP. The engineering of TTP mutants with modulating binding properties can have potential impact at industrial level for easier purification of single tocopherols from vitamin E mixtures coming from natural oils or synthetic processes. Moreover, the identification of a -tocopherol selective TTP offers the possibility to challenge the hypotheses for the evolutionary development of a mechanism for -tocopherol selection in omnivorous animals.     

Molecular characterization of the substance P*neurokinin-1 receptor complex: development of an experimentally based model

Molecular models for the interaction of substance P (SP) with its G protein-coupled receptor, the neurokinin-1 receptor (NK-1R), have been developed. The ligand.receptor complex is based on experimental data from a series of photoaffinity labeling experiments and spectroscopic structural studies of extracellular domains of the NK-1R. Using the ligand/receptor contact points derived from incorporation of photolabile probes (p-benzoylphenylalanine (Bpa)) into SP at positions 3, 4, and 8 and molecular dynamics simulations, the topological arrangement of SP within the NK-1R is explored. The model incorporates the structural features, determined by high resolution NMR studies, of the second extracellular loop (EC2), containing contact points Met(174) and Met(181), providing important experimentally based conformational preferences for the simulations. Extensive molecular dynamics simulations were carried out to probe the nature of the two contact points identified for the Bpa(3)SP analogue (Bremer, A. A., Leeman, S. E., and Boyd, N. D. (2001) J. Biol. Chem. 276, 22857-22861), examining modes of ligand binding in which the contact points are fulfilled sequentially or simultaneously. The resulting ligand.receptor complex has the N terminus of SP projecting toward transmembrane helix (TM) 1 and TM2, exposed to the solvent. The C terminus of SP is located in proximity to TM5 and TM6, deeper into the central core of the receptor. The central portion of the ligand, adopting a helical loop conformation, is found to align with the helices of the central regions EC2 and EC3, forming important interactions with both of these extracellular domains. The model developed here allows for atomic insight into the biochemical data currently available and guides targeting of future experiments to probe specific ligand/receptor interactions and thereby furthers our understanding of the functioning of this important neuropeptide system.     

New insights into the GABA(A) receptor structure and orthosteric ligand binding: receptor modeling guided by experimental data

GABA(A) receptors (GABA(A)Rs) are ligand gated chloride ion channels that mediate overall inhibitory signaling in the CNS. A detailed understanding of their structure is important to gain insights in, e.g., ligand binding and functional properties of this pharmaceutically important target. Homology modeling is a necessary tool in this regard because experimentally determined structures are lacking. Here we present an exhaustive approach for creating a high quality model of the α(1)β(2)γ(2) subtype of the GABA(A)R ligand binding domain, and we demonstrate its usefulness in understanding details of orthosteric ligand binding. The model was constructed by using multiple templates and by incorporation of knowledge from biochemical/pharmacological experiments. It was validated on the basis of objective energy functions, its ability to account for available residue specific information, and its stability in molecular dynamics (MD) compared with that of the two homologous crystal structures. We then combined the model with extensive structure-activity relationships available from two homologous series of orthosteric GABA(A)R antagonists to create a detailed hypothesis for their binding modes. Excellent agreement with key experimental data was found, including the ability of the model to accommodate and explain a previously developed pharmacophore model. A coupling to agonist binding was thereby established and discussed in relation to activation mechanisms. Our results highlight the importance of critical evaluation and optimization of each step in the homology modeling process. The approach taken here can greatly aid in increasing the understanding of GABA(A)Rs and related receptors where structural insight is limited and reliable models are difficult to obtain.     

Access Path to the Ligand Binding Pocket May Play a Role in Xenobiotics Selection by AhR

Understanding of multidrug binding at the atomic level would facilitate drug design and strategies to modulate drug metabolism, including drug transport, oxidation, and conjugation. Therefore we explored the mechanism of promiscuous binding of small molecules by studying the ligand binding domain, the PAS-B domain of the aryl hydrocarbon receptor (AhR). Because of the low sequence identities of PAS domains to be used for homology modeling, structural features of the widely employed HIF-2α and a more recent suitable template, CLOCK were compared. These structures were used to build AhR PAS-B homology models. We performed molecular dynamics simulations to characterize dynamic properties of the PAS-B domain and the generated conformational ensembles were employed in in silico docking. In order to understand structural and ligand binding features we compared the stability and dynamics of the promiscuous AhR PAS-B to other PAS domains exhibiting specific interactions or no ligand binding function. Our exhaustive in silico binding studies, in which we dock a wide spectrum of ligand molecules to the conformational ensembles, suggest that ligand specificity and selection may be determined not only by the PAS-B domain itself, but also by other parts of AhR and its protein interacting partners. We propose that ligand binding pocket and access channels leading to the pocket play equally important roles in discrimination of endogenous molecules and xenobiotics.     

Molecular dynamics simulations of apo and holo forms of fatty acid binding protein 5 and cellular retinoic acid binding protein II reveal highly mobile protein,retinoic acid ligand,and water molecules

Structural and dynamic properties from a series of 300 ns molecular dynamics, MD, simulations of two intracellular lipid binding proteins, iLBPs, (Fatty Acid Binding Protein 5, FABP5, and Cellular Retinoic Acid Binding Protein II, CRABP-II) in both the apo form and when bound with retinoic acid reveal a high degree of protein and ligand flexibility. The ratio of FABP5 to CRABP-II in a cell may determine whether it undergoes natural apoptosis or unrestricted cell growth in the presence of retinoic acid. As a result, FABP5 is a promising target for cancer therapy. The MD simulations presented here reveal distinct differences in the two proteins and provide insight into the binding mechanism. CRABP-II is a much larger, more flexible protein that closes upon ligand binding, where FABP5 transitions to an open state in the holo form. The traditional understanding obtained from crystal structures of the gap between two β-sheets of the β-barrel common to iLBPs and the α-helix cap that forms the portal to the binding pocket is insufficient for describing protein conformation (open vs. closed) or ligand entry and exit. When the high degree of mobility between multiple conformations of both the ligand and protein are examined via MD simulation, a new mode of ligand motion that improves understanding of binding dynamics is revealed.     

Unmasking ligand binding motifs: identification of a chemokine receptor motif by NMR studies of antagonist peptides

Determining the critical structural features a ligand must possess in order to bind to its receptor is of key importance to the understanding of vital biological processes and to the rational design of small molecule therapeutics to modulate receptor function. We have developed a general strategy for determining such ligand binding motifs using low temperature NMR structures of peptides with the desired receptor binding properties. This approach has been successfully applied to determine a binding motif for the chemokine receptor CXCR4. The motif identified provides a detailed guide for the design of small molecule antagonists against CXCR4, which are much sought after to aid in the treatment of a number of conditions including human immunodeficiency virus type 1 infection and a variety of cancers.     

Molecular mechanism of substrate specificity in the bacterial neutral amino acid transporter LeuT

The recently published X-ray structure of LeuT, a Na(+)/Cl(-)-dependent neurotransmitter transporter, has provided fresh impetus to efforts directed at understanding the molecular principles governing specific neurotransmitter transport. The combination of the LeuT crystal structure with the results of molecular simulations enables the functional data on specific binding and transport to be related to molecular structure. All-atom FEP and molecular dynamics (MD) simulations of LeuT embedded in an explicit membrane were performed alongside a decomposition analysis to dissect the molecular determinants of the substrate specificity of LeuT. It was found that the ligand must be in a zwitterionic (ZW) form to bind tightly to the transporter. The theoretical results on the absolute binding-free energies for leucine, alanine, and glycine show that alanine can be a potent substrate for LeuT, although leucine is preferred, which is consistent with the recent experimental data (Singh et al., Nature 2007;448:952-956). Furthermore, LeuT displays robust specificity for leucine over glycine. Interestingly, the ability of LeuT to discriminate between substrates relies on the dynamics of residues that form its binding pocket (e.g., F253 and Q250) and the charged side chains (R30-D404) from a second coordination shell. The water-mediated R30-D404 salt bridge is thought to be part of the extracellular (EC) gate of LeuT. The introduction of a polar ligand such as glycine to the water-depleted binding pocket of LeuT gives rise to structural rearrangements of the R30-D404-Q250 hydrogen-bonding network and leads to increased hydration of the binding pocket. Conformational changes associated with the broken hydrogen bond between Q250 and R30 are shown to be important for tight and selective ligand binding to LeuT.     

Retinal Conformation Changes Rhodopsin’s Dynamic Ensemble

G protein-coupled receptors are vital membrane proteins that allosterically transduce biomolecular signals across the cell membrane. However, the process by which ligand binding induces protein conformation changes is not well understood biophysically. Rhodopsin, the mammalian dim-light receptor, is a unique test case for understanding these processes because of its switch-like activity; the ligand, retinal, is bound throughout the activation cycle, switching from inverse agonist to agonist after absorbing a photon. By contrast, the ligand-free opsin is outside the activation cycle and may behave differently. We find that retinal influences rhodopsin dynamics using an ensemble of all-atom molecular dynamics simulations that in aggregate contain 100 μs of sampling. Active retinal destabilizes the inactive state of the receptor, whereas the active ensemble was more structurally homogenous. By contrast, simulations of an active-like receptor without retinal present were much more heterogeneous than those containing retinal. These results suggest allosteric processes are more complicated than a ligand inducing protein conformational changes or simply capturing a shifted ensemble as outlined in classic models of allostery.     

Molecular Determinants of Epidermal Growth Factor Binding: A Molecular Dynamics Study

The epidermal growth factor receptor (EGFR) is a member of the receptor tyrosine kinase family that plays a role in multiple cellular processes. Activation of EGFR requires binding of a ligand on the extracellular domain to promote conformational changes leading to dimerization and transphosphorylation of intracellular kinase domains. Seven ligands are known to bind EGFR with affinities ranging from sub-nanomolar to near micromolar dissociation constants. In the case of EGFR, distinct conformational states assumed upon binding a ligand is thought to be a determining factor in activation of a downstream signaling network. Previous biochemical studies suggest the existence of both low affinity and high affinity EGFR ligands. While these studies have identified functional effects of ligand binding, high-resolution structural data are lacking. To gain a better understanding of the molecular basis of EGFR binding affinities, we docked each EGFR ligand to the putative active state extracellular domain dimer and 25.0 ns molecular dynamics simulations were performed. MM-PBSA/GBSA are efficient computational approaches to approximate free energies of protein-protein interactions and decompose the free energy at the amino acid level. We applied these methods to the last 6.0 ns of each ligand-receptor simulation. MM-PBSA calculations were able to successfully rank all seven of the EGFR ligands based on the two affinity classes: EGF>HB-EGF>TGF-α>BTC>EPR>EPG>AR. Results from energy decomposition identified several interactions that are common among binding ligands. These findings reveal that while several residues are conserved among the EGFR ligand family, no single set of residues determines the affinity class. Instead we found heterogeneous sets of interactions that were driven primarily by electrostatic and Van der Waals forces. These results not only illustrate the complexity of EGFR dynamics but also pave the way for structure-based design of therapeutics targeting EGF ligands or the receptor itself.     

DynaPred: a structure and sequence based method for the prediction of MHC class I binding peptide sequences and conformations

MOTIVATION: The binding of endogenous antigenic peptides to MHC class I molecules is an important step during the immunologic response of a host against a pathogen. Thus, various sequence- and structure-based prediction methods have been proposed for this purpose. The sequence-based methods are computationally efficient, but are hampered by the need of sufficient experimental data and do not provide a structural interpretation of their results. The structural methods are data-independent, but are quite time-consuming and thus not suited for screening of whole genomes. Here, we present a new method, which performs sequence-based prediction by incorporating information obtained from molecular modeling. This allows us to perform large databases screening and to provide structural information of the results. RESULTS: We developed a SVM-trained, quantitative matrix-based method for the prediction of MHC class I binding peptides, in which the features of the scoring matrix are energy terms retrieved from molecular dynamics simulations. At the same time we used the equilibrated structures obtained from the same simulations in a simple and efficient docking procedure. Our method consists of two steps: First, we predict potential binders from sequence data alone and second, we construct protein-peptide complexes for the predicted binders. So far, we tested our approach on the HLA-A0201 allele. We constructed two prediction models, using local, position-dependent (DynaPred(POS)) and global, position-independent (DynaPred) features. The former model outperformed the two sequence-based methods used in our evaluation; the latter shows a much higher generalizability towards other alleles than the position-dependent models. The constructed peptide structures can be refined within seconds to structures with an average backbone RMSD of 1.53 A from the corresponding experimental structures.     

Cell,Isoform, and Environment Factors Shape Gradients and Modulate Chemotaxis

Chemokine gradient formation requires multiple processes that include ligand secretion and diffusion, receptor binding and internalization, and immobilization of ligand to surfaces. To understand how these events dynamically shape gradients and influence ensuing cell chemotaxis, we built a multi-scale hybrid agent-based model linking gradient formation, cell responses, and receptor-level information. The CXCL12/CXCR4/CXCR7 signaling axis is highly implicated in metastasis of many cancers. We model CXCL12 gradient formation as it is impacted by CXCR4 and CXCR7, with particular focus on the three most highly expressed isoforms of CXCL12. We trained and validated our model using data from an in vitro microfluidic source-sink device. Our simulations demonstrate how isoform differences on the molecular level affect gradient formation and cell responses. We determine that ligand properties specific to CXCL12 isoforms (binding to the migration surface and to CXCR4) significantly impact migration and explain differences in in vitro chemotaxis data. We extend our model to analyze CXCL12 gradient formation in a tumor environment and find that short distance, steep gradients characteristic of the CXCL12-γ isoform are effective at driving chemotaxis. We highlight the importance of CXCL12-γ in cancer cell migration: its high effective affinity for both extracellular surface sites and CXCR4 strongly promote CXCR4+ cell migration. CXCL12-γ is also more difficult to inhibit, and we predict that co-inhibition of CXCR4 and CXCR7 is necessary to effectively hinder CXCL12-γ-induced migration. These findings support the growing importance of understanding differences in protein isoforms, and in particular their implications for cancer treatment.     

A hydrophobic loop in acyl-CoA binding protein is functionally important for binding to palmitoyl-coenzyme A: a molecular dynamics study

Acyl-CoA binding protein (ACBP) plays a key role in lipid metabolism, interacting via a partly unknown mechanism with high affinity with long chain fatty acyl-CoAs (LCFA-CoAs). At present there is no study of the microscopic way ligand binding is accomplished. We analyzed this process by molecular dynamics (MDs) simulations. We proposed a computational model of ligand, able to reproduce some evidence from nuclear magnetic resonance (NMR) data, quantitative time resolved fluorometry and X-ray crystallography. We found that a hydrophobic loop, not in the active site, is important for function. Besides, multiple sequence alignment shows hydrophobicity (and not the residues itselves) conservation.     

Quaternary solution structures of galectins-1, -3, and -7

Galectins are a growing family of animal lectins with functions in growth regulation and cell adhesion that bind beta-Gal residues in oligosaccharides. Evidence indicates that some of the biological properties of galectins are due to their cross-linking activities with multivalent glycoconjugate receptors. Therefore determination of the quaternary solution structures of these proteins is important in understanding their structure-function properties. The present study reports analytical sedimentation velocity and equilibrium data for galectins-1, -3, and -7 in the absence and presence of bound LacNAc, the natural ligand epitope. Galectin-1 from bovine heart and recombinant human galectin-7 were found to be stable dimers by both methods. In contrast, recombinant murine galectin-3, as well as its proteolytical derived C-terminal domain, are predominantly monomeric. The presence of LacNAc at concentrations sufficient to fully saturate the proteins had no significant effect on either the weight average molecular weight determined by sedimentation equilibrium or the hydrodynamic properties determined from sedimentation velocity experiments. These results show that binding of a monovalent ligand does not affect oligomerization of these galectins.     

A universal thermodynamic approach to analyze biomolecular binding experiments

Binding processes of any kind can be characterized as an association of a given ligand with some binding factor. This includes macromolecules as well as supramolecular aggregates such as micelles or membranes. The underlying molecular binding mechanism may be more or less complicated due to various intermediate steps (involving for instance conformational changes, aggregation, cooperativity, etc.). A sensible discussion of possible binding models naturally calls for a model-independent access to basic thermodynamic properties. The present contribution will demonstrate how this can quite generally be accomplished by a pertinent processing of properly selected experimental data. The method requires a series of titration measurements comprising the use of variable amounts of both the ligand and the binding factor. It leads to a linear mass conservation plot (i.e. amount of the ligand vs. a matching amount of the binding factor) whose slope and ordinate intercept are equal to the binding ratio (i.e. bound ligand per binding factor) and the free ligand concentration, respectively. This establishes the specific binding isotherm. The approach also reveals latent structurally determined features of the applied physical measuring signal. A number of examples including specific binding, unspecific adsorption and insertion in two-dimensional molecular films will illustrate the methodology.