Microsecond-timescale Simulations Research Articles

Caging Polycations: Effect of Increasing Confinement on the Modes of Interaction of Spermidine3+ With DNA Double Helices.

Polyamines have important roles in the modulation of the cellular function and are ubiquitous in cells. The polyamines putrescine2+, spermidine3+, and spermine4+ represent the most abundant organic counterions of the negatively charged DNA in the cellular nucleus. These polyamines are known to stabilize the DNA structure and, depending on their concentration and additional salt composition, to induce DNA aggregation, which is often referred to as condensation. However, the modes of interactions of these elongated polycations with DNA and how they promote condensation are still not clear. In the present work, atomistic molecular dynamics (MD) computer simulations of two DNA fragments surrounded by spermidine3+ (Spd3+) cations were performed to study the structuring of Spd3+ “caged” between DNA molecules. Microsecond time scale simulations, in which the parallel DNA fragments were constrained at three different separations, but allowed to rotate axially and move naturally, provided information on the conformations and relative orientations of surrounding Spm3+ cations as a function of DNA-DNA separation. Novel geometric criteria allowed for the classification of DNA-Spd3+ interaction modes, with special attention given to Spd3+ conformational changes in the space between the two DNA molecules (caged Spd3+). This work shows how changes in the accessible space, or confinement, around DNA affect DNA-Spd3+ interactions, information fundamental to understanding the interactions between DNA and its counterions in environments where DNA is compacted, e.g. in the cellular nucleus.

Mechanism of RNA recognition by a Musashi RNA-binding protein.

The Musashi RNA-binding proteins (RBPs) regulate translation of target mRNAs and maintenance of cell stemness and tumorigenesis. Musashi-1 (MSI1), long considered as an intestinal and neural stem cell marker, has been more recently found to be over expressed in many cancers. It has served as an important drug target for treating acute myeloid leukemia and solid tumors such as ovarian, colorectal and bladder cancer. One of the reported binding targets of MSI1 is Numb, a negative regulator of the Notch signaling. However, the dynamic mechanism of Numb RNA binding to MSI1 remains unknown, largely hindering effective drug design targeting this critical interaction. Here, we have performed extensive all-atom microsecond-timescale simulations using a robust Gaussian accelerated molecular dynamics (GaMD) method, which successfully captured multiple times of spontaneous and highly accurate binding of the Numb RNA from bulk solvent to the MSI1 protein target site. GaMD simulations revealed that Numb RNA binding to MSI1 involved largely induced fit in both the RNA and protein. The simulations also identified important low-energy intermediate conformational states during RNA binding, in which Numb interacted mainly with the β2-β3 loop and C terminus of MSI1. The mechanistic understanding of RNA binding obtained from our GaMD simulations is expected to facilitate rational structure-based drug design targeting MSI1 and other RBPs.

Microscopic Description of Protein-RNA Interactions in Nucleoprotein Condensates

All eukaryotic cells use membrane-less compartments to spatially confine critical cellular processes, thereby increasing their throughput. A growing number of studies suggest that such compartmentalization occurs in response to cellular signals via a physical mechanism akin to liquid-liquid phase separation. Fused in sarcoma (FUS) is one of the best-studied proteins known to undergo a liquid-liquid phase separation, forming a biological condensate that recruits specific biomolecules implicated in DNA repair. Recent experiments have found the physical properties of such condensates to depend sensitively on the presence of nucleic acids in the condensate and on point mutations to the FUS sequence. A number of such point mutations have been implicated in neurodegenerative diseases including amyotrophic lateral sclerosis and dementia, for which there is no cure. By carrying out long time scale, all-atom molecular dynamics simulations of FUS protein systems, this study provides the first atomic-level insights into the structure of the biological condensate. The study's primary objective was to identify specific interactions within the FUS proteins, as well as between the proteins and RNA, that control the condensate's viscosity. Microsecond time scale simulations of a single FUS protein in solution and of an aggregate of two FUS proteins revealed the formation of persistent arginine-tyrosine contacts that were previously reported to orchestrate FUS aggregation. Uracil RNA homopolymers of varying lengths were then introduced within the FUS aggregate and simulated using the all-atom molecular dynamics method. The simulations have shown that binding of FUS to RNA depends on the length of the RNA, suggesting a pathway for RNA import into the condensate. Finally, we have applied the steered molecular dynamics protocol to pull bound RNA molecules out of the condensate, which has elucidated the critical molecular interactions holding the condensates together.

Microsecond-timescale simulations suggest 5-HT–mediated preactivation of the 5-HT3A serotonin receptor

Aided by efforts to improve their speed and efficiency, molecular dynamics (MD) simulations provide an increasingly powerful tool to study the structure-function relationship of pentameric ligand-gated ion channels (pLGICs). However, accurate reporting of the channel state and observation of allosteric regulation by agonist binding with MD remains difficult due to the timescales necessary to equilibrate pLGICs from their artificial and crystalized conformation to a more native, membrane-bound conformation in silico. Here, we perform multiple all-atom MD simulations of the homomeric 5-hydroxytryptamine 3A (5-HT3A) serotonin receptor for 15 to 20 μs to demonstrate that such timescales are critical to observe the equilibration of a pLGIC from its crystalized conformation to a membrane-bound conformation. These timescales, which are an order of magnitude longer than any previous simulation of 5-HT3A, allow us to observe the dynamic binding and unbinding of 5-hydroxytryptamine (5-HT) (i.e., serotonin) to the binding pocket located on the extracellular domain (ECD) and allosteric regulation of the transmembrane domain (TMD) from synergistic 5-HT binding. While these timescales are not long enough to observe complete activation of 5-HT3A, the allosteric regulation of ion gating elements by 5-HT binding is indicative of a preactive state, which provides insight into molecular mechanisms that regulate channel activation from a resting state. This mechanistic insight, enabled by microsecond-timescale MD simulations, will allow a careful examination of the regulation of pLGICs at a molecular level, expanding our understanding of their function and elucidating key structural motifs that can be targeted for therapeutic regulation.

Leveraging local MP2 to reduce basis set superposition errors: An efficient first-principles based force-field for carbon dioxide.

Pairwise additive model potentials for CO2 were developed by fitting to gradients computed with the local second order Møller Plesset Perturbation theory (LMP2) method, with and without consideration of 3-body dispersion using adaptive force matching. Without fitting to experiments, all models gave good predictions of properties of CO2, such as the density-temperature diagram, diffusion constants, and radial distribution functions. For the prediction of vibrational spectra, the inclusion of a bond-bond coupling term has been shown to be important. The CO2 models developed only have pairwise additive terms, thus allowing microsecond time scale simulations to be performed with practical computational cost. LMP2 performed significantly better than second order Møller Plesset Perturbation theory (MP2) for the development of the CO2 model. This is attributed to the appreciable reduction in the basis set superposition error when the localized method was used. It is argued that LMP2 is a more appropriate method than MP2 for force matching for systems where the basis set superposition error is large.

Characterization of Zebrafish Cardiac and Slow Skeletal Troponin C Paralogs by MD Simulation and ITC

Zebrafish, as a model for teleost fish, have two paralogous troponin C (TnC) genes that are expressed in the heart differentially in response to temperature acclimation. Upon Ca2+ binding, TnC changes conformation and exposes a hydrophobic patch that interacts with troponin I and initiates cardiac muscle contraction. Teleost-specific TnC paralogs have not yet been functionally characterized. In this study we have modeled the structures of the paralogs using molecular dynamics simulations at 18°C and 28°C and calculated the different Ca2+-binding properties between the teleost cardiac (cTnC or TnC1a) and slow-skeletal (ssTnC or TnC1b) paralogs through potential-of-mean-force calculations. These values are compared with thermodynamic binding properties obtained through isothermal titration calorimetry (ITC). The modeled structures of each of the paralogs are similar at each temperature, with the exception of helix C, which flanks the Ca2+ binding site; this region is also home to paralog-specific sequence substitutions that we predict have an influence on protein function. The short timescale of the potential-of-mean-force calculation precludes the inclusion of the conformational change on the ΔG of Ca2+ interaction, whereas the ITC analysis includes the Ca2+ binding and conformational change of the TnC molecule. ITC analysis has revealed that ssTnC has higher Ca2+ affinity than cTnC for Ca2+ overall, whereas each of the paralogs has increased affinity at 28°C compared to 18°C. Microsecond-timescale simulations have calculated that the cTnC paralog transitions from the closed to the open state more readily than the ssTnC paralog, an unfavorable transition that would decrease the ITC-derived Ca2+ affinity while simultaneously increasing the Ca2+ sensitivity of the myofilament. We propose that the preferential expression of cTnC at lower temperatures increases myofilament Ca2+ sensitivity by this mechanism, despite the lower Ca2+ affinity that we have measured by ITC.

Probing the Antimicrobial Action of Polymyxin B1 and Melittin VIA Coarse-Grained Molecular Dynamics Simulations

Many pathogenic Gram-negative bacteria are progressively acquiring resistance to the available spectrum of antibiotic drugs. Cationic antimicrobial peptides, which are known to exhibit rapid in vitro bactericidal activity, and a low propensity for resistance development, have been identified as promising candidates for the development of novel antibacterial agents. Rational development of therapeutics based on antimicrobial agents requires a molecular-level understanding of their mechanism of action. To contribute to this endeavour, we have studied two very different antimicrobial peptides, the cyclic lipopeptide, polymyxin B1 and the single helix peptide, melittin. We have performed a series of coarse-grained molecular dynamics simulations to explore how both peptides interact with the outer membrane of Gram-negative bacteria. The microsecond timescale simulations reveal very different patterns of behaviour of the two peptides and also show the very marked effect of having an asymmetric membrane. The latter observation is particularly crucial as this key feature of the in vivo bacterial outer membrane is often neglected in computational and experimental studies.

The Structural Basis for Activation and Inhibition of ZAP-70 Kinase Domain.

ZAP–70 (Zeta-chain-associated protein kinase 70) is a tyrosine kinase that interacts directly with the activated T-cell receptor to transduce downstream signals, and is hence a major player in the regulation of the adaptive immune response. Dysfunction of ZAP–70 causes selective T cell deficiency that in turn results in persistent infections. ZAP–70 is activated by a variety of signals including phosphorylation of the kinase domain (KD), and binding of its regulatory tandem Src homology 2 (SH2) domains to the T cell receptor. The present study investigates molecular mechanisms of activation and inhibition of ZAP–70 via atomically detailed molecular dynamics simulation approaches. We report microsecond timescale simulations of five distinct states of the ZAP–70 KD, comprising apo, inhibited and three phosphorylated variants. Extensive analysis of local flexibility and correlated motions reveal crucial transitions between the states, thus elucidating crucial steps in the activation mechanism of the ZAP–70 KD. Furthermore, we rationalize previously observed staurosporine-bound crystal structures, suggesting that whilst the KD superficially resembles an “active-like” conformation, the inhibitor modulates the underlying protein dynamics and restricts it in a compact, rigid state inaccessible to ligands or cofactors. Finally, our analysis reveals a novel, potentially druggable pocket in close proximity to the activation loop of the kinase, and we subsequently use its structure in fragment-based virtual screening to develop a pharmacophore model. The pocket is distinct from classical type I or type II kinase pockets, and its discovery offers promise in future design of specific kinase inhibitors, whilst mutations in residues associated with this pocket are implicated in immunodeficiency in humans.

Accurately modeling nanosecond protein dynamics requires at least microseconds of simulation.

Advances in hardware and algorithms have greatly extended the timescales accessible to molecular simulation. This article assesses whether such long timescale simulations improve our ability to calculate order parameters that describe conformational heterogeneity on ps-ns timescales or if force fields are now a limiting factor. Order parameters from experiment are compared with order parameters calculated in three different ways from simulations ranging from 10 ns to 100 μs in length. Importantly, bootstrapping is employed to assess the variability in results for each simulation length. The results of 10-100 ns timescale simulations are highly variable, possibly explaining the variation in levels of agreement between simulation and experiment in published works examining different proteins. Fortunately, microsecond timescale simulations improve both the accuracy and precision of calculated order parameters, at least so long as the full exponential fit or truncated average approximation is used instead of the common long-time limit approximation. The improved precision of these long timescale simulations allows a statistically sound comparison of a number of modern force fields, such as Amber03, Amber99sb-ILDN, and Charmm27. While there is some variation between these force fields, they generally give very similar results for sufficiently long simulations. The fact that so much simulation is required to precisely capture ps-ns timescale processes suggests that extremely extensive simulations are required for slower processes. Advanced sampling techniques could aid greatly, however, such methods will need to maintain accurate kinetics if they are to be of value for calculating dynamical properties like order parameters. © 2015 Wiley Periodicals, Inc.

Small Fdu Strand Exhibits Salt-Dependent Stability

As shown in our previous work, F10 (a 10mer of Fdu - 5-fluoro-2′-deoxyuridine-5′-O-monophosphate) has a 40% likelihood of folding into a stable hairpin structure when Magnesium ions are used to neutralize the simulated system. Investigations of other counter ions revealed that F10 folds less than 10% of the time when neutralized with ions with a +1 charge and 40% of the time when neutralized with ions with a +2 charge. When simulated with biologically relevant combinations of K1+, Na1+ and Mg2+ ions, F10 folds into this hairpin structure approximately 10% of the time; however, once the strand of Fdu enters this folded state, it is 96% likely to stay in this kinetic trap. Here we present the data from microsecond timescale simulations and hypotheses on F10's behavior in the presence of various counter ions.

Molecular Dynamics Study of Ion Conduction and Selectivity in a Prokaryotic Ion Channel

Since the publication of crystal structure of NavAb, the first of a (prokaryotic) voltage gated sodium channel, several computational studies have been aimed at determining the mechanisms of ion selectivity and ion conduction in voltage gated sodium channels. We provide a two-part study, well-converged free energy surfaces involving one, two and three ions provide results consistent with microsecond timescale simulations of ion conduction over concentration and voltages. The position of ions in the pore and cavity correlate to coordination number and side chain orientation, showing, as suggested by Chakrabarti et al. (PNAS (2013) 110, 11331). More surprising, presence of an ion in the aqueous cavity beneath the selectivity filter was sufficient to influence glutamate conformation.Because the crystal structure of NavAb was in a closed pore conformation, we truncated the S5 and S6 helices and restrained these outer helices harmonically in a membrane represented by supporting lattice of neon. The turrets, pore helices and selectivity filter were free to move. Conductances were in the range of the experimental values; however, our studies show no preference for sodium conduction over potassium conduction at physiological conditions (−200mV, 0.15 M salt). While sodium conductance varied little with respect to concentration and presence of potassium in solution, stronger potassium selectivity is observed at 1 M and in mixed solutions. Our results are consistent with the observation of anomalous mole fraction effect in NaChBac by Finol-Urdaneta et al. (JGP (2014) 143, 157-171), though not sufficient to confirm such an effect. Additionally, we provide mechanisms for ion conduction showing a greater variety of states and transitions occupied during potassium conduction, and note that at higher concentrations of salt the mechanism of conduction is not as clearly cut.

Role of Desolvation in Thermodynamics and Kinetics of Ligand Binding to a Kinase.

Computer simulations are used to determine the free energy landscape for the binding of the anticancer drug Dasatinib to its src kinase receptor and show that before settling into a free energy basin the ligand must surmount a free energy barrier. An analysis based on using both the ligand-pocket separation and the pocket-water occupancy as reaction coordinates shows that the free energy barrier is a result of the free energy cost for almost complete desolvation of the binding pocket. The simulations further show that the barrier is not a result of the reorganization free energy of the binding pocket. Although a continuum solvent model gives the location of free energy minima, it is not able to reproduce the intermediate free energy barrier. Finally, it is shown that a kinetic model for the on rate constant in which the ligand diffuses up to a doorway state and then surmounts the desolvation free energy barrier is consistent with published microsecond time-scale simulations of the ligand binding kinetics for this system [Shaw, D. E. et al. J. Am. Chem. Soc.2011, 133, 9181−918321545110].

Ff14ipq: A Self-Consistent Force Field for Condensed-Phase Simulations of Proteins.

We present the ff14ipq force field, implementing the previously published IPolQ charge set for simulations of complete proteins. Minor modifications to the charge derivation scheme and van der Waals interactions between polar atoms are introduced. Torsion parameters are developed through a generational learning approach, based on gas-phase MP2/cc-pVTZ single-point energies computed of structures optimized by the force field itself rather than the quantum benchmark. In this manner, we sacrifice information about the true quantum minima in order to ensure that the force field maintains optimal agreement with the MP2/cc-pVTZ benchmark for the ensembles it will actually produce in simulations. A means of making the gas-phase torsion parameters compatible with solution-phase IPolQ charges is presented. The ff14ipq model is an alternative to ff99SB and other Amber force fields for protein simulations in programs that accommodate pair-specific Lennard–Jones combining rules. The force field gives strong performance on α-helical and β-sheet oligopeptides as well as globular proteins over microsecond time scale simulations, although it has not yet been tested in conjunction with lipid and nucleic acid models. We show how our choices in parameter development influence the resulting force field and how other choices that may have appeared reasonable would actually have led to poorer results. The tools we developed may also aid in the development of future fixed-charge and even polarizable biomolecular force fields.

Event detection and sub‐state discovery from biomolecular simulations using higher‐order statistics: Application to enzyme adenylate kinase

Biomolecular simulations at millisecond and longer time-scales can provide vital insights into functional mechanisms. Because post-simulation analyses of such large trajectory datasets can be a limiting factor in obtaining biological insights, there is an emerging need to identify key dynamical events and relating these events to the biological function online, that is, as simulations are progressing. Recently, we have introduced a novel computational technique, quasi-anharmonic analysis (QAA) (Ramanathan et al., PLoS One 2011;6:e15827), for partitioning the conformational landscape into a hierarchy of functionally relevant sub-states. The unique capabilities of QAA are enabled by exploiting anharmonicity in the form of fourth-order statistics for characterizing atomic fluctuations. In this article, we extend QAA for analyzing long time-scale simulations online. In particular, we present HOST4MD--a higher-order statistical toolbox for molecular dynamics simulations, which (1) identifies key dynamical events as simulations are in progress, (2) explores potential sub-states, and (3) identifies conformational transitions that enable the protein to access those sub-states. We demonstrate HOST4MD on microsecond timescale simulations of the enzyme adenylate kinase in its apo state. HOST4MD identifies several conformational events in these simulations, revealing how the intrinsic coupling between the three subdomains (LID, CORE, and NMP) changes during the simulations. Further, it also identifies an inherent asymmetry in the opening/closing of the two binding sites. We anticipate that HOST4MD will provide a powerful and extensible framework for detecting biophysically relevant conformational coordinates from long time-scale simulations.

Agonist Dynamics and Conformational Selection during Microsecond Simulations of the A2A Adenosine Receptor

The G-protein-coupled receptors (GPCRs) are a ubiquitous family of signaling proteins of exceptional pharmacological importance. The recent publication of structures of several GPCRs cocrystallized with ligands of differing activity offers a unique opportunity to gain insight into their function. To that end, we performed microsecond-timescale simulations of the A2A adenosine receptor bound to either of two agonists, adenosine or UK432097. Our data suggest that adenosine is highly dynamic when bound to A2A, in stark contrast to the case with UK432097. Remarkably, adenosine finds an alternate binding pose in which the ligand is inverted relative to the crystal structure, forming relatively stable interactions with helices I and II. Our observations suggest new experimental tests to validate our predictions and deepen our understanding of GPCR signaling. Overall, our data suggest an intriguing hypothesis: that the 100- to 1000-fold greater efficacy of UK432097 relative to adenosine arises because UK432097 stabilizes a much tighter neighborhood of active conformations, which manifests as a greater likelihood of G-protein activation per unit time.

Ligand Dynamics during Microsecond Simulations of the A2a Adenosine Receptor

Thanks to their diversity of functions and widespread expression in all tissues, activation of G-protein-coupled receptors (GPCRs) upon agonist binding plays a critical role in cell signaling in many contexts. The GPCRs are therefore one of the most important classes of protein as drug targets. Recently, the crystal structure of the A2a adenosine receptor bound to agonist adenosine (G. Lebon et al., Nature 474, 521 (2011)) or UK432097 (F. Xu et al., Science 332, 322 (2011)) has been reported, offering insight into receptor A2a ligand binding and activation. Using the crystal structures as a starting point, we performed microsecond-timescale simulations of the A2a receptors bound to adenosine or UK432097, as well as the antagonist ZM241385. Our data suggest that adenosine is highly dynamic when bound to A2a, which contrasts starkly with the larger ligands UK432097 and ZM241385. Analyzing the simulation data in the context of published mutagensis experimental data suggests that the extra molecular scaffold of the larger ligands serves to stabilize key interactions, perhaps rationalizing the relative efficacy of different agonists.

Insights into the folding pathway of the Engrailed Homeodomain protein using replica exchange molecular dynamics simulations

Protein folding studies were carried out by performing microsecond time scale simulations on the ultrafast/fast folding protein Engrailed Homeodomain (EnHD) from Drosophila melanogaster. It is a three-helix bundle protein consisting of 54 residues (PDB ID: 1ENH). The positions of the helices are 8–20 (Helix I), 26–36 (Helix II) and 40–53 (Helix III). The second and third helices together form a Helix-Turn-Helix (HTH) motif which belongs to the family of DNA binding proteins. The molecular dynamics (MD) simulations were performed using replica exchange molecular dynamics (REMD). REMD is a method that involves simulating a protein at different temperatures and performing exchanges at regular time intervals. These exchanges were accepted or rejected based on the Metropolis criterion. REMD was performed using the AMBER FF03 force field with the generalised Born solvation model for the temperature range 286–373 K involving 30 replicas. The extended conformation of the protein was used as the starting structure. A simulation of 600 ns per replica was performed resulting in an overall simulation time of 18 μs. The protein was seen to fold close to the native state with backbone root mean square deviation (RMSD) of 3.16 Å. In this low RMSD structure, the Helix I was partially formed with a backbone RMSD of 3.37 Å while HTH motif had an RMSD of 1.81 Å. Analysis suggests that EnHD folds to its native structure via an intermediate in which the HTH motif is formed. The secondary structure development occurs first followed by tertiary packing. The results were in good agreement with the experimental findings.

Automated Event Detection and Activity Monitoring in Long Molecular Dynamics Simulations.

Events of scientific interest in molecular dynamics (MD) simulations, including conformational changes, folding transitions, and translocations of ligands and reaction products, often correspond to high-level structural rearrangements that alter contacts between molecules or among different parts of a molecule. Due to advances in computer architecture and software, MD trajectories representing such structure-changing events have become easier to generate, but the length of these trajectories poses a challenge to scientific interpretation and analysis. In this paper, we present automated methods for the detection of potentially important structure-changing events in long MD trajectories. In contrast with traditional tools for the analysis of such trajectories, our methods provide a detailed report of broken and formed contacts that aids in the identification of specific time-dependent side-chain interactions. Our approach employs a coarse-grained representation of amino acid side chains, a contact metric based on higher order generalizations of Delaunay tetrahedralization, techniques for detecting significant shifts in the resulting contact time series, and a new kernel-based measure of contact alteration activity. The analysis methods we describe are incorporated in a newly developed package, called TimeScapes, which is freely available and compatible with trajectories generated by a variety of popular MD programs. Tests based on actual microsecond time scale simulations demonstrate that the package can be used to efficiently detect and characterize important conformational changes in realistic protein systems.

Automated Event Detection and Activity Monitoring in Long Time-Scale Molecular Dynamics

Molecular dynamics trajectories of biological systems contain many events of great scientific interest, including conformational transitions, folding processes, and translocations of ligands and reaction products. In proteins these events often correspond to high-level tertiary or quaternary structure rearrangements, which alter the contacts between amino acid residues. Due to advances in computer architecture and software, molecular dynamics trajectories representing such structure-changing events have become easy to generate, but their length complicates scientific interpretation. The goal of this work is to simplify an important part of the analysis workflow (and to complement traditional visual inspection) by automating the mining of long trajectories. We present several new algorithms and implementation techniques that enable the detection of significant structure-changing events in a molecular dynamics trajectory. These algorithms include a coarse graining of side chain contacts, a contact metric based on higher-order generalizations of the Delaunay tetrahedralization, and median filters for detecting significant shifts in the ensemble mean of the resulting time series. We have also developed numerical techniques for suppressing trivial re-crossing events and a new kernel-based estimator of the contact alteration activity. These methods will be disseminated in a newly developed package, “TimeScapes,” which is compatible with molecular dynamics trajectories generated from any of a variety of popular simulation programs. Tests on microsecond time scale simulations suggest that the implementation is efficient and requires very little parameterization. The analysis provides a detailed listing of broken and formed contacts, and reliably detects allosteric and folding transitions, as well as stable intermediates, in the protein dynamics.