Search

Intrinsically disordered proteins (IDPs) are implicated in many human diseases. They have generally not been amenable to conventional structure-based drug design, however, because their intrinsic conformational variability has precluded an atomic-level understanding of their binding to small molecules. Here we present long-timescale, atomic-level molecular dynamics (MD) simulations of monomeric α-synuclein (an IDP whose aggregation is associated with Parkinson’s disease) binding the small-molecule drug fasudil in which the observed protein-ligand interactions were found to be in good agreement with previously reported NMR chemical shift data. In our simulations, fasudil, when bound, favored certain charge-charge and π-stacking interactions near the C terminus of α-synuclein, but tended not to form these interactions simultaneously, rather breaking one of these interactions and forming another nearby (a mechanism we term dynamic shuttling). Further simulations with small molecules chosen to modify these interactions yielded binding affinities and key structural features of binding consistent with subsequent NMR experiments, suggesting the potential for MD-based strategies to facilitate the rational design of small molecules that bind with disordered proteins.

Significance Most proteins associate with other proteins to function, forming complexes that are central to almost all physiological processes. Determining the structures of these complexes and understanding how they associate are problems of fundamental importance. Using long-timescale molecular dynamics simulations, some performed using a new enhanced sampling method, we observed spontaneous association and dissociation of five protein–protein systems to and from their experimentally determined native complexes. By analyzing the simulations of these five systems, which include members of diverse structural and functional classes, we are able to draw general mechanistic conclusions about protein association., Despite the biological importance of protein–protein complexes, determining their structures and association mechanisms remains an outstanding challenge. Here, we report the results of atomic-level simulations in which we observed five protein–protein pairs repeatedly associate to, and dissociate from, their experimentally determined native complexes using a molecular dynamics (MD)–based sampling approach that does not make use of any prior structural information about the complexes. To study association mechanisms, we performed additional, conventional MD simulations, in which we observed numerous spontaneous association events. A shared feature of native association for these five structurally and functionally diverse protein systems was that if the proteins made contact far from the native interface, the native state was reached by dissociation and eventual reassociation near the native interface, rather than by extensive interfacial exploration while the proteins remained in contact. At the transition state (the conformational ensemble from which association to the native complex and dissociation are equally likely), the protein–protein interfaces were still highly hydrated, and no more than 20% of native contacts had formed.

A quantitative characterization of the binding properties of drug fragments to a target protein is an important component of a fragment-based drug discovery program. Fragments typically have a weak binding affinity, however, making it challenging to experimentally characterize key binding properties, including binding sites, poses, and affinities. Direct simulation of the binding equilibrium by molecular dynamics (MD) simulations can provide a computational route to characterize fragment binding, but this approach is so computationally intensive that it has thus far remained relatively unexplored. Here, we perform MD simulations of sufficient length to observe several different fragments spontaneously and repeatedly bind to and unbind from the protein FKBP, allowing the binding affinities, on- and off-rates, and relative occupancies of alternative binding sites and alternative poses within each binding site to be estimated, thereby illustrating the potential of long time scale MD as a quantitative tool for fragment-based drug discovery. The data from the long time scale fragment binding simulations reported here also provide a useful benchmark for testing alternative computational methods aimed at characterizing fragment binding properties. As an example, we calculated binding affinities for the same fragments using a standard free energy perturbation approach and found that the values agreed with those obtained from the fragment binding simulations within statistical error.

Significance Bruton’s tyrosine kinase (Btk) activation on the cell membrane is critical for B cell proliferation and development, and Btk inhibition is a promising treatment for several hematologic cancers and autoimmune diseases. Here, we examine Btk activation using the results of long-timescale molecular dynamics simulations. In our simulations, Btk lipid-binding modules dimerized on the membrane in a single predominant conformation. We observed that the phospholipid PIP3—in addition to its expected role of recruiting Btk to the membrane—allosterically mediated dimer formation and stability by binding at two novel sites. Our results provide strong evidence that PIP3-mediated dimerization of Btk at the cell membrane is a critical step in Btk activation and suggest a potential approach to allosteric Btk inhibitor development., Bruton’s tyrosine kinase (Btk) is critical for B cell proliferation and activation, and the development of Btk inhibitors is a vigorously pursued strategy for the treatment of various B cell malignancies. A detailed mechanistic understanding of Btk activation has, however, been lacking. Here, inspired by a previous suggestion that Btk activation might depend on dimerization of its lipid-binding PH–TH module on the cell membrane, we performed long-timescale molecular dynamics simulations of membrane-bound PH–TH modules and observed that they dimerized into a single predominant conformation. We found that the phospholipid PIP3 stabilized the dimer allosterically by binding at multiple sites, and that the effects of PH–TH mutations on dimer stability were consistent with their known effects on Btk activity. Taken together, our simulation results strongly suggest that PIP3-mediated dimerization of Btk at the cell membrane is a critical step in Btk activation.

Molecular Basis of Ligand Dissociation from the Adenosine A2A ReceptorDong Guo, Albert C. Pan, Ron O. Dror, Tamara Mocking, Rongfang Liu, Laura H. Heitman, David E. Shaw and Adriaan P. IJzermanMolecular Pharmacology May 2016, 89 (5) 485-491; DOI: https://doi.org/10.1124/mol.115.102657AbstractHow drugs dissociate from their targets is largely unknown. We investigated the molecular basis of this process in the adenosine A2A receptor (A2AR), a prototypical G protein–coupled receptor (GPCR). Through kinetic radioligand binding experiments, we characterized mutant receptors selected based on molecular dynamic simulations of the antagonist ZM241385 dissociating from the A2AR. We discovered mutations that dramatically altered the ligand’s dissociation rate despite only marginally influencing its binding affinity, demonstrating that even receptor features with little contribution to affinity may prove critical to the dissociation process. Our results also suggest that ZM241385 follows a multistep dissociation pathway, consecutively interacting with distinct receptor regions, a mechanism that may also be common to many other GPCRs.

The pathways that G protein-coupled receptor (GPCR) ligands follow as they bind to or dissociate from their receptors are largely unknown. Protease-activated receptor-1 (PAR1) is a GPCR activated by intramolecular binding of a tethered agonist peptide that is exposed by thrombin cleavage. By contrast, the PAR1 antagonist vorapaxar is a lipophilic drug that binds in a pocket almost entirely occluded from the extracellular solvent. The binding and dissociation pathway of vorapaxar is unknown. Starting with the crystal structure of vorapaxar bound to PAR1, we performed temperature-accelerated molecular dynamics simulations of ligand dissociation. In the majority of simulations, vorapaxar exited the receptor laterally into the lipid bilayer through openings in the transmembrane helix (TM) bundle. Prior to full dissociation, vorapaxar paused in metastable intermediates stabilized by interactions with the receptor and lipid headgroups. Derivatives of vorapaxar with alkyl chains predicted to extend between TM6 and TM7 into the lipid bilayer inhibited PAR1 with apparent on rates similar to that of the parent compound in cell signaling assays. These data are consistent with vorapaxar binding to PAR1 via a pathway that passes between TM6 and TM7 from the lipid bilayer, in agreement with the most consistent pathway observed by molecular dynamics. While there is some evidence of entry of the ligand into rhodopsin and lipid-activated GPCRs from the cell membrane, our study provides the first such evidence for a peptide-activated GPCR and suggests that metastable intermediates along drug binding and dissociation pathways can be stabilized by specific interactions between lipids and the ligand.

Bruton’s tyrosine kinase (Btk) is critical for B-cell proliferation and activation, and the development of Btk inhibitors is a vigorously pursued strategy for the treatment of various B-cell malignancies. A detailed mechanistic understanding of Btk activation has, however, been lacking. Here, inspired by a previous suggestion that Btk activation might depend on dimerization of its lipid-binding PH-TH module on the cell membrane, we performed long-timescale molecular dynamics simulations of membrane-bound PH-TH modules and observed that they dimerized into a single predominant conformation. We found that the phospholipid PIP3 stabilized the dimer allosterically by binding at multiple sites, and that the effects of PH-TH mutations on dimer stability were consistent with their known effects on Btk activity. Taken together, our simulation results strongly suggest that PIP3-mediated dimerization of Btk at the cell membrane is a critical step in Btk activation.

Binding modes and molecular mechanisms of several allosteric modulators of a prototypical G-protein-coupled receptor are revealed using atomic-level simulations and validated by the rational design of a modulator with substantially altered effects. A third of clinically used drugs elicit their biological effects via a G-protein-coupled receptor (GPCR), usually by binding at the orthosteric (primary ligand-binding) site, in competition with the ligands that naturally regulate receptor signalling. The design of small molecules able to selectively modulate a GPCR by binding to an allosteric site is a desirable goal, but difficult to achieve because neither the binding modes nor the molecular mechanisms of such molecules are known. In this manuscript, the authors used molecular dynamics simulations, with experimental validation, to determine where and how structurally diverse allosteric modulators bind to the M2 muscarinic acetylcholine receptor, which is essential for the physiological control of cardiovascular function. Despite substantial structural diversity of the small molecule ligands, the molecules all formed cation-π interactions with clusters of aromatic residues in the receptor's extracellular vestibule, about 15 A from the ligand-binding pocket. These findings may facilitate the rational design of allosteric modulators targeting muscarinic and related GPCRs. The design of G-protein-coupled receptor (GPCR) allosteric modulators, an active area of modern pharmaceutical research, has proved challenging because neither the binding modes nor the molecular mechanisms of such drugs are known1,2. Here we determine binding sites, bound conformations and specific drug–receptor interactions for several allosteric modulators of the M2 muscarinic acetylcholine receptor (M2 receptor), a prototypical family A GPCR, using atomic-level simulations in which the modulators spontaneously associate with the receptor. Despite substantial structural diversity, all modulators form cation–π interactions with clusters of aromatic residues in the receptor extracellular vestibule, approximately 15 A from the classical, ‘orthosteric’ ligand-binding site. We validate the observed modulator binding modes through radioligand binding experiments on receptor mutants designed, on the basis of our simulations, either to increase or to decrease modulator affinity. Simulations also revealed mechanisms that contribute to positive and negative allosteric modulation of classical ligand binding, including coupled conformational changes of the two binding sites and electrostatic interactions between ligands in these sites. These observations enabled the design of chemical modifications that substantially alter a modulator’s allosteric effects. Our findings thus provide a structural basis for the rational design of allosteric modulators targeting muscarinic and possibly other GPCRs.

Application of a specific stimulus opens the intracellular gate of a K(+) channel (activation), yielding a transient period of ion conduction until the selectivity filter spontaneously undergoes a conformational change towards a non-conductive state (inactivation). Removal of the stimulus closes the gate and allows the selectivity filter to interconvert back to its conductive conformation (recovery). Given that the structural differences between the conductive and inactivated filter are very small, it is unclear why the recovery process can take up to several seconds. The bacterial K(+) channel KcsA from Streptomyces lividans can be used to help elucidate questions about channel inactivation and recovery at the atomic level. Although KcsA contains only a pore domain, without voltage-sensing machinery, it has the structural elements necessary for ion conduction, activation and inactivation. Here we reveal, by means of a series of long molecular dynamics simulations, how the selectivity filter is sterically locked in the inactive conformation by buried water molecules bound behind the selectivity filter. Potential of mean force calculations show how the recovery process is affected by the buried water molecules and the rebinding of an external K(+) ion. A kinetic model deduced from the simulations shows how releasing the buried water molecules can stretch the timescale of recovery to seconds. This leads to the prediction that reducing the occupancy of the buried water molecules by imposing a high osmotic stress should accelerate the rate of recovery, which was verified experimentally by measuring the recovery rate in the presence of a 2-molar sucrose concentration.

G-protein-coupled receptors (GPCRs) can modulate diverse signaling pathways, often in a ligand-specific manner. The full range of functionally relevant GPCR conformations is poorly understood. Here, we use NMR spectroscopy to characterize the conformational dynamics of the transmembrane core of the β(2)-adrenergic receptor (β(2)AR), a prototypical GPCR. We labeled β(2)AR with (13)CH(3)ε-methionine and obtained HSQC spectra of unliganded receptor as well as receptor bound to an inverse agonist, an agonist, and a G-protein-mimetic nanobody. These studies provide evidence for conformational states not observed in crystal structures, as well as substantial conformational heterogeneity in agonist- and inverse-agonist-bound preparations. They also show that for β(2)AR, unlike rhodopsin, an agonist alone does not stabilize a fully active conformation, suggesting that the conformational link between the agonist-binding pocket and the G-protein-coupling surface is not rigid. The observed heterogeneity may be important for β(2)AR's ability to engage multiple signaling and regulatory proteins.

Acetylcholine, the first neurotransmitter to be identified, exerts many of its physiological actions via activation of a family of G-protein-coupled receptors (GPCRs) known as muscarinic acetylcholine receptors (mAChRs). Although the five mAChR subtypes (M1-M5) share a high degree of sequence homology, they show pronounced differences in G-protein coupling preference and the physiological responses they mediate. Unfortunately, despite decades of effort, no therapeutic agents endowed with clear mAChR subtype selectivity have been developed to exploit these differences. We describe here the structure of the G(q/11)-coupled M3 mAChR ('M3 receptor', from rat) bound to the bronchodilator drug tiotropium and identify the binding mode for this clinically important drug. This structure, together with that of the G(i/o)-coupled M2 receptor, offers possibilities for the design of mAChR subtype-selective ligands. Importantly, the M3 receptor structure allows a structural comparison between two members of a mammalian GPCR subfamily displaying different G-protein coupling selectivities. Furthermore, molecular dynamics simulations suggest that tiotropium binds transiently to an allosteric site en route to the binding pocket of both receptors. These simulations offer a structural view of an allosteric binding mode for an orthosteric GPCR ligand and provide additional opportunities for the design of ligands with different affinities or binding kinetics for different mAChR subtypes. Our findings not only offer insights into the structure and function of one of the most important GPCR families, but may also facilitate the design of improved therapeutics targeting these critical receptors.

A third of marketed drugs act by binding to a G-protein-coupled receptor (GPCR) and either triggering or preventing receptor activation. Although recent crystal structures have provided snapshots of both active and inactive functional states of GPCRs, these structures do not reveal the mechanism by which GPCRs transition between these states. Here we propose an activation mechanism for the β 2 -adrenergic receptor, a prototypical GPCR, based on atomic-level simulations in which an agonist-bound receptor transitions spontaneously from the active to the inactive crystallographically observed conformation. A loosely coupled allosteric network, comprising three regions that can each switch individually between multiple distinct conformations, links small perturbations at the extracellular drug-binding site to large conformational changes at the intracellular G-protein-binding site. Our simulations also exhibit an intermediate that may represent a receptor conformation to which a G protein binds during activation, and suggest that the first structural changes during receptor activation often take place on the intracellular side of the receptor, far from the drug-binding site. By capturing this fundamental signaling process in atomic detail, our results may provide a foundation for the design of drugs that control receptor signaling more precisely by stabilizing specific receptor conformations.

n/a

The coupled interplay between activation and inactivation gating is a functional hallmark of K(+) channels. This coupling has been experimentally demonstrated through ion interaction effects and cysteine accessibility, and is associated with a well defined boundary of energetically coupled residues. The structure of the K(+) channel KcsA in its fully open conformation, in addition to four other partial channel openings, richly illustrates the structural basis of activation-inactivation gating. Here, we identify the mechanistic principles by which movements on the inner bundle gate trigger conformational changes at the selectivity filter, leading to the non-conductive C-type inactivated state. Analysis of a series of KcsA open structures suggests that, as a consequence of the hinge-bending and rotation of the TM2 helix, the aromatic ring of Phe 103 tilts towards residues Thr 74 and Thr 75 in the pore-helix and towards Ile 100 in the neighbouring subunit. This allows the network of hydrogen bonds among residues Trp 67, Glu 71 and Asp 80 to destabilize the selectivity filter, allowing entry to its non-conductive conformation. Mutations at position 103 have a size-dependent effect on gating kinetics: small side-chain substitutions F103A and F103C severely impair inactivation kinetics, whereas larger side chains such as F103W have more subtle effects. This suggests that the allosteric coupling between the inner helical bundle and the selectivity filter might rely on straightforward mechanical deformation propagated through a network of steric contacts. Average interactions calculated from molecular dynamics simulations show favourable open-state interaction-energies between Phe 103 and the surrounding residues. We probed similar interactions in the Shaker K(+) channel where inactivation was impaired in the mutant I470A. We propose that side-chain rearrangements at position 103 mechanically couple activation and inactivation in KcsA and a variety of other K(+) channels.

An approach to find transition pathways in complex systems is presented. The method, which is related to the string method in collective variables of Maragliano et al. [J.Chem. Phys. 125:024106 (2006)], is conceptually simple and straightforward to implement. It consists in refining a putative transition path in the multi-dimensional space supported by a set of collective variables using the average dynamic drift of those variables. This drift is estimated on-the-fly via swarms of short unbiased trajectories started at different points along the path. Successive iterations of this algorithm, which can be naturally distributed over many computer nodes with negligible inter-processor communication, refine an initial trial path toward the most probable transition path (MPTP) between two stable basins. The method is first tested by determining the pathway for the C7eq to C7ax transition in an all-atom model of the alanine dipeptide in vacuum, which has been studied previously with the string method in collective variables. A transition path is found with a committor distribution peaked at 1/2 near the free energy maximum, in accord with previous results. Lastly, the method is applied to the allosteric conformational change in the nitrogen regulatory protein C (NtrC), represented here with a two-state elastic network model. Even though more than 550 collective variables are used to describe the conformational change, the path converges rapidly. Again, the committor distribution is found to be peaked around 1/2 near the free energy maximum between the two stable states, confirming that a genuine transition state has been localized in this complex multi-dimensional system.

Structurally elucidating transition pathways between protein conformations gives deep mechanistic insight into protein behavior but is typically difficult. Unbiased molecular dynamics (MD) simulations provide one solution, but their computational expense is often prohibitive, motivating the development of enhanced sampling methods that accelerate conformational changes in a given direction, embodied in a collective variable. The accuracy of such methods is unclear for complex protein transitions, because obtaining unbiased MD data for comparison is difficult. Here, we use long-time scale, unbiased MD simulations of epidermal growth factor receptor kinase deactivation as a complex biological test case for two widely used methods-steered molecular dynamics (SMD) and the string method. We found that common collective variable choices, based on the root-mean-square deviation (RMSD) of the entire protein, prevented the methods from producing accurate paths, even in SMD simulations on the time scale of the unbiased transition. Using collective variables based on the RMSD of the region of the protein known to be important for the conformational change, however, enabled both methods to provide a more accurate description of the pathway in a fraction of the simulation time required to observe the unbiased transition.

We have used a combination of neutron scattering experiments and Monte Carlo simulations to study the initial stages of first-order phase transitions. We focus on quenches wherein the nascent phase is formed by homogeneous nucleation, and we approach the spinodal, i.e., the quench depth at which the original phase becomes unstable. In this regime, we show how critical nuclei sizes are determined from neutron scattering structure factors. Prevailing thought is that the size of the critical nucleus should increase with increasing quench depth and diverge at the spinodal. To the contrary, our experiments and simulations indicate that the critical nucleus size decreases monotonically as quench depth is increased and is finite at the spinodal.

Despite the biological importance of protein-protein complexes, determining their structures and association mechanisms remains an outstanding challenge. Here, we report the results of atomic-level simulations in which we observed five protein-protein pairs repeatedly associate to, and dissociate from, their experimentally determined native complexes using a new molecular dynamics (MD)-based sampling approach that does not make use of any prior structural information about the complexes. To study association mechanisms, we performed additional, conventional MD simulations, in which we observed numerous spontaneous association events. A shared feature of native association for these five structurally and functionally diverse protein systems was that if the proteins made contact far from the native interface, the native state was reached by dissociation and eventual re-association near the native interface, rather than by extensive interfacial exploration while the proteins remained in contact. At the transition state (the conformational ensemble from which association to the native complex and dissociation are equally likely), the protein-protein interfaces were still highly hydrated, and no more than 20% of native contacts had formed.

All-atom molecular dynamics (MD) simulation is a valuable technique for providing detailed information about the dynamics of biomolecules, but its computational expense can often be prohibitive. This limitation has motivated the development of “enhanced sampling” simulation methods - purely algorithmic changes to conventional MD that aim to accelerate the sampling of configurational states. Although many such methods are available, relatively few systematic studies of their performance have been done. Quantitative claims about the performance of enhanced sampling simulations are typically limited to (i) comparisons with conventional MD simulations of small, model systems, which may present qualitatively different sampling challenges than complex biological systems, or (ii) comparisons with experimental data, which may be complicated by discrepancies between the actual experimental conditions and the modeling of those experimental conditions in the enhanced sampling simulation, including errors in the physical model, or force field, used in simulation. An effective alternative approach to quantifying performance of enhanced sampling methods would be to compare the results they obtain on complex biological systems directly to those obtained by conventional MD simulations using the same force field. Here, we use long-timescale conventional MD simulations to assess the performance of certain commonly used enhanced sampling methods in accelerating the sampling of protein conformational changes, protein folding, and ligand binding.

Understanding the nature of the glass transition--the dramatic slowing of dynamics and eventual emergence of a disordered solid from a cooling liquid--is a fundamental challenge in physical science. A central characteristic of glass-forming liquids is a non-exponential main relaxation process. The extent of deviation from exponential relaxation typically becomes more pronounced on cooling. Theories that predict a growth of spatially heterogeneous dynamics as temperature is lowered can explain these observations. In apparent contradiction to these theories, however, some experiments suggest that certain substances--notably including the intensely studied molecular glass-former ortho-terphenyl (OTP)--have a main relaxation process whose shape is essentially temperature independent, even though other observables predicted to be correlated with the degree of dynamical heterogeneity are temperature dependent. Here we report the first simulations based on an atomistic model of OTP that reach equilibrium at temperatures well into the supercooled regime. We first show that the results of these simulations are in reasonable quantitative agreement with experimental data for several basic properties over a wide range of temperatures. We then focus on rotational relaxation, finding nearly exponential behavior at high temperatures with clearly increasing deviations as temperature is lowered. The much weaker temperature dependence observed in light-scattering experiments also emerges from the same simulation data when we calculate correlation functions similar to those probed experimentally; this highlights the diversity of temperature dependencies that can be obtained with different probes. Further analysis suggests that the temperature insensitivity observed in the light-scattering experiments stems from the dependence of these measurements on internal as well as rotational molecular motion. Within the temperature range of our OTP simulations, our results strongly suggest that this archetypal glass-former behaves as anticipated by theories of the glass transition that predict increasing non-exponentiality with cooling, and our simulations thus strengthen the evidence supporting such theories.

Special-purpose computing hardware can provide significantly better performance and power efficiency for certain applications than general-purpose processors. Even within a single application area, however, a special-purpose machine can be far more valuable if it is capable of efficiently supporting a number of different computational methods that, taken together, expand the machine's functionality and range of applicability. We have previously described a massively parallel special-purpose supercomputer, called Anton, and have shown that it executes traditional molecular dynamics simulations orders of magnitude faster than the previous state of the art. Here, we describe how we extended Anton's software to support a more diverse set of methods, allowing scientists to simulate a broader class of biological phenomena at extremely high speeds. Key elements of our approach, which exploits Anton's tightly integrated hardwired pipelines and programmable cores, are applicable to the hardware and software design of various other specialized or heterogeneous parallel computing platforms.

It is increasingly appreciated that the rates at which drugs associate with and dissociate from receptors--the binding kinetics--directly impact drug efficacy and safety. The molecular determinants of drug-receptor binding kinetics remain poorly understood, however, especially when compared with the well-known factors that affect binding affinity. The rational modulation of kinetics during lead optimization thus remains challenging. We review some of the key factors thought to control drug-receptor binding kinetics at the molecular level--molecular size, conformational fluctuations, electrostatic interactions and hydrophobic effects--and discuss several possible approaches for the rational design of drugs with desired binding kinetics.

SummaryIn K+ channels, rearrangements of the pore outer vestibule have been associated with C-type inactivation gating. Paradoxically, the crystal structure of Open/C-type inactivated KcsA suggests these movements to be modest in magnitude. In this study, we show that under physiological conditions, the KcsA outer vestibule undergoes relatively large dynamic rearrangements upon inactivation. External Cd2+ enhances the rate of C-type inactivation in an cysteine mutant (Y82C) via metal-bridge formation. This effect is not present in a non-inactivating mutant (E71A/Y82C). Tandem dimer and tandem tetramer constructs of equivalent cysteine mutants in KcsA and Shaker K+ channels demonstrate that these Cd2+ metal bridges are formed only between adjacent subunits. This is well supported by molecular dynamics simulations. Based on the crystal structure of Cd2+-bound Y82C-KcsA in the closed state, together with electron paramagnetic resonance distance measurements in the KcsA outer vestibule, we suggest that subunits must dynamically come in close proximity as the channels undergo inactivation.

How drugs bind to their receptors—from initial association, through drug entry into the binding pocket, to adoption of the final bound conformation, or “pose”—has remained unknown, even for G-protein-coupled receptor modulators, which constitute one-third of all marketed drugs. We captured this pharmaceutically critical process in atomic detail using the first unbiased molecular dynamics simulations in which drug molecules spontaneously associate with G-protein-coupled receptors to achieve final poses matching those determined crystallographically. We found that several beta blockers and a beta agonist all traverse the same well-defined, dominant pathway as they bind to the β 1 - and β 2 -adrenergic receptors, initially making contact with a vestibule on each receptor’s extracellular surface. Surprisingly, association with this vestibule, at a distance of 15 Å from the binding pocket, often presents the largest energetic barrier to binding, despite the fact that subsequent entry into the binding pocket requires the receptor to deform and the drug to squeeze through a narrow passage. The early barrier appears to reflect the substantial dehydration that takes place as the drug associates with the vestibule. Our atomic-level description of the binding process suggests opportunities for allosteric modulation and provides a structural foundation for future optimization of drug–receptor binding and unbinding rates.

The amount of ionic current flowing through K+ channels is determined by the interplay between two separate time-dependent processes: activation and inactivation gating. Activation is concerned with the stimulus-dependent opening of the main intracellular gate, whereas inactivation is a spontaneous conformational transition of the selectivity filter toward a nonconductive state occurring on a variety of timescales. A recent analysis of multiple x-ray structures of open and partially open KcsA channels revealed the mechanism by which movements of the inner activation gate, formed by the inner helices from the four subunits of the pore domain, bias the conformational changes at the selectivity filter toward a nonconductive inactivated state. This analysis highlighted the important role of Phe103, a residue located along the inner helix, near the hinge position associated with the opening of the intracellular gate. In the present study, we use free energy perturbation molecular dynamics simulations (FEP/MD) to quantitatively elucidate the thermodynamic basis for the coupling between the intracellular gate and the selectivity filter. The results of the FEP/MD calculations are in good agreement with experiments, and further analysis of the repulsive, van der Waals dispersive, and electrostatic free energy contributions reveals that the energetic basis underlying the absence of inactivation in the F103A mutation in KcsA is the absence of the unfavorable steric interaction occurring with the large Ile100 side chain in a neighboring subunit when the intracellular gate is open and the selectivity filter is in a conductive conformation. Macroscopic current analysis shows that the I100A mutant indeed relieves inactivation in KcsA, but to a lesser extent than the F103A mutant.

Modal-gating shifts represent an effective regulatory mechanism by which ion channels control the extent and time course of ionic fluxes. Under steady-state conditions, the K+ channel KcsA displays three distinct gating modes, high-Po, low-Po and a high-frequency flicker mode, each with about an order of magnitude difference in their mean open times. In KcsA, the hydrogen bond network between Glu71, Asp80 and Trp67 that surrounds the selectivity filter has been shown to regulate C-type inactivation, with the Glu71-Asp80 pair having the strongest influence on selectivity filter stability. Here, we show that in the absence of C-type inactivation, mutations at the pore-helix position Glu71 unmask a series of kinetically distinct modes of gating in a side-chain-specific way which mirror those seen in wild-type channels. Results from high-resolution crystal structures along with molecular dynamic simulations suggest that specific interactions in the side-chain network surrounding the selectivity filter, in concert with ion occupancy, alter the relative stability of pre-existing conformational states of the pore. These findings highlight the key role of the selectivity filter in regulating modal gating behavior in K+ channels.

An efficient method is proposed for building Markov models with discrete states able to accurately describe the slow relaxation of a complex system with two stable conformations. First, the reaction pathway described by a set of collective variables between the two stable states is determined using the string method with swarms of trajectories. Then, short trajectories are initiated at different points along this pathway to build the state-to-state transition probability matrix. It is shown, using a model system, how this strategy makes it possible to use trajectories that are significantly shorter than the slowest relaxation time to efficiently build a reliable and accurate Markov model. Extensions of the method to multiple pathways, as well as some common pitfalls arising from poorly relaxed paths or an inappropriate choice of collective variables, are illustrated and discussed.

Dynamical four-point susceptibilities measure the extent of spatial correlations in the dynamics of glass forming systems. We show how these susceptibilities depend on the length scales that necessarily form part of their definition. The behaviour of these susceptibilities is estimated by means of an analysis in terms of renewal processes within the context of dynamic facilitation. The analytic results are confirmed by numerical simulations of an atomistic model glass-former, and of two kinetically constrained models. Hence we argue that the scenario predicted by the dynamic facilitation approach is generic., 10 pages, 8 figures, v3 has minor revisions

We study dynamical heterogeneity and glassy dynamics in a kinetically constrained lattice gas model which has both translational and rotational degrees of freedom. We find that the rotational diffusion constant tracks the structural relaxation time as density is increased whereas the translational diffusion constant exhibits a strong decoupling. We investigate distributions of exchange and persistence times for both the rotational and translational degrees of freedom and compare our results on the distributions of rotational exchange times to recent single molecule studies., 7 pages, 5 figures

We study dynamical heterogeneity and growing dynamical lengthscales in two kinetically constrained models, namely, the one- and two-vacancy assisted triangular lattice gases. One of the models is a strong glassformer and the other is a fragile glassformer. Both exhibit heterogeneous dynamics with broadly distributed timescales as seen in the distribution of persistence times. We show that the Stokes-Einstein relation is violated, to a greater degree in the fragile glassformer, and show how this violation is related to dynamic heterogeneity. We extract dynamical lengthscales from structure factors of mobile particles and show, quantitatively, the growth of this lengthscale as density increases. We comment on how the scaling of lengths and times in these models relates to that in facilitated spin models of glasses., 9 pages, 13 figures, updated fig. 12 and added fig. 13

We present an interpolated kinetically constrained lattice gas model which exhibits a transition from fragile to strong supercooled liquid behavior. We find non-monotonic decoupling that is due to this crossover and is seen in experiment., 3 pages, 3 figures, in press at Chem. Phys. Chem

Reactive pathways to nucleation in a three-dimensional Ising model at 60% of the critical temperature are studied using transition path sampling of single spin flip Monte Carlo dynamics. Analysis of the transition state ensemble (TSE) indicates that the critical nuclei are rough and anisotropic. The TSE, projected onto the free energy surface characterized by cluster size, N, and surface area, S, indicates the significance of other variables in addition to these two traditional reaction coordinates for nucleation. The transmission coefficient along N is ~ 0.35, and this reduction of the transmission coefficient from unity is explained in terms of the stochastic nature of the dynamic model., In press at the Journal of Physical Chemistry B, 7 pages, 8 figures

The design of G protein-coupled receptor (GPCR) allosteric modulators, an active area of modern pharmaceutical research, has proven challenging because neither the binding modes nor the molecular mechanisms of such drugs are known. Here we determined binding sites, bound conformations, and specific drug-receptor interactions for several allosteric modulators of the M2 muscarinic acetylcholine receptor, a prototypical Family A GPCR, using atomic-level simulations in which the modulators spontaneously associated with the receptor (Nature, in press). Despite substantial structural diversity, all modulators formed cation-π interactions with clusters of aromatic residues in the receptor extracellular vestibule, ∼15 A from the classical, “orthosteric” ligand-binding site. We validated the observed modulator binding modes through radioligand binding experiments on receptor mutants designed, on the basis of our simulations, to either increase or decrease modulator affinity. Simulations also revealed mechanisms that contribute to positive and negative allosteric modulation of classical ligand binding, including coupled conformational changes of the two binding sites and electrostatic interactions between ligands in these sites. These observations enabled the design of chemical modifications that significantly altered a modulator's allosteric effects. Our findings thus provide a structural basis for the rational design of allosteric modulators targeting muscarinic and possibly other GPCRs.

Potassium channels control the flow of ions across cell membranes with gating mechanisms involving conformational changes at the intracellular gate and the selectivity filter. Opening of the intracellular gate via external stimuli (activation) results in a transient period of conduction before the selectivity filter undergoes a conformational change, which constricts the permeation pathway (inactivation). When the applied stimulus is removed and the lower gate closes (deactivation), the filter slowly resets to a conductive conformation (recovery from inactivation). Using the KcsA channel as a prototypical model system to examine these issues, a combination of computer simulation (all-atom free energy and potential of mean force computations as well as transition pathway determination using the string method with swarms-of-trajectories) and experiment (electrophysiology and X-ray crystallography) is used to provide new insight into the microscopic mechanism that underlies inactivation and recovery from inactivation. An ion binding event is revealed as a crucial step in resetting the inactive filter during the recovery from inactivation.

Understanding the mechanism of conformational changes in macromolecules requires the knowledge of the intermediate states. A version of the string method, which uses multiple short dynamics trajectories to propagate the pathway, was recently developed by Pan et al. Here we use data from swarms of trajectories calculated at discrete points in phase space to interpolate the average displacement and variance at arbitrary points. This is tested on model potentials using statistics from actual swarms of trajectories. We use the interpolated parameters to compute the Markovian propagators from one point on the transition path to the next. We use them to obtain a time-dependent action of a path, which can be optimized to produce the highest probability pathway. We describe the optimization protocol and demonstrate that in artificial flat potentials the existing string method cannot correct problems such as loops in the initial path, while the new method produces the correct pathway (Figure shows pathway in 2D potential). We further illustrate the utility of our method by applying it to protein conformational transitions, such as the KcsA potassium channel, and comparing its performance to existing transition pathway methods.View Large Image | View Hi-Res Image | Download PowerPoint Slide

Intrinsically disordered proteins (IDPs) are implicated in many human diseases. They have generally not been amenable to conventional structure-based drug design, however, because their intrinsic conformational variability has precluded an atomic-level understanding of their binding to small molecules. Here we present long-time-scale, atomic-level molecular dynamics (MD) simulations of monomeric α-synuclein (an IDP whose aggregation is associated with Parkinson's disease) binding the small-molecule drug fasudil in which the observed protein-ligand interactions were found to be in good agreement with previously reported NMR chemical shift data. In our simulations, fasudil, when bound, favored certain charge-charge and π-stacking interactions near the C terminus of α-synuclein but tended not to form these interactions simultaneously, rather breaking one of these interactions and forming another nearby (a mechanism we term dynamic shuttling). Further simulations with small molecules chosen to modify these interactions yielded binding affinities and key structural features of binding consistent with subsequent NMR experiments, suggesting the potential for MD-based strategies to facilitate the rational design of small molecules that bind with disordered proteins.

We have known the structure for the closed-state of a potassium channel pore domain (PD) for more than a decade. However, major progress in understanding the molecular basis for activation and inactivation gating in K-channels had to wait until high-resolution structural information of the channel in the open state became available. Recently, we solved the structure KcsA in it fully open conformation, as well four others partial openings, which richly illustrated the channel activation-inactivation pathway. Analysis of these open structures suggested that residue F103 in TM2 interacts with the c-terminal end of the pore helix, compressing the pitch of its first helical turn. As a consequence, the distance between E71-D80 side chains is shortened, strengthening the carboxyl-carboxylate interaction that leads to a non-conductive conformation of the selectivity filter. Perturbation mutagenesis at position 103, affected gating kinetics as predicted from our structural analysis: small side chain substitutions F103A and F103C severely impaired inactivation kinetics, suggesting an allosteric coupling between the inner helical bundle and the selectivity filter. Free energy calculations show strong open state interaction-energies between F103 and surrounding residues. Similar interactions were probed in the Shaker K-channel by mutating highly conserved I470, equivalent to F103, to a smaller side chain. In the mutant I470A, inactivation was abrogated, suggesting that a similar mechanism underlies inactivation coupling in eukaryotic potassium channels. A crystallography study of these mutants in the open KcsA will be reported.

The potassium ion channel, KcsA, gates the passage of ions through cell membranes in response to a change in pH. Recent experimental results have demonstrated the existence of two gates in KcsA: an intracellular gate and a gate at the selectivity filter. Lowering the pH opens the intracellular gate allowing ions to pass. After a period of time, however, the channel inactivates by constricting the selectivity filter and impeding the flow of ions even though the bottom gate remains open. We have used path-based molecular dynamics simulations to probe the detailed mechanism of this phenomenon by finding dynamical pathways to inactivation in KcsA. We have computed free energies and rates of inactivation that agree with recent experimental results. We also provide a molecular rationalization for the coupling between the opening and closing of the lower gate and the inactivation of the selectivity filter.