Multi-microsecond Molecular Dynamics Simulations Research Articles

Excitatory and inhibitory D-serine binding to the NMDA receptor.

N-methyl-D-aspartate receptors (NMDARs) uniquely require binding of two different neurotransmitter agonists for synaptic transmission. D-serine and glycine bind to one subunit, GluN1, while glutamate binds to the other, GluN2. These agonists bind to the receptor's bi-lobed ligand-binding domains (LBDs), which close around the agonist during receptor activation. To better understand the unexplored mechanisms by which D-serine contributes to receptor activation, we performed multi-microsecond molecular dynamics simulations of the GluN1/GluN2A LBD dimer with free D-serine and glutamate agonists. Surprisingly, we observed D-serine binding to both GluN1 and GluN2A LBDs, suggesting that D-serine competes with glutamate for binding to GluN2A. This mechanism is confirmed by our electrophysiology experiments, which show that D-serine is indeed inhibitory at high concentrations. Although free energy calculations indicate that D-serine stabilizes the closed GluN2A LBD, its inhibitory behavior suggests that it either does not remain bound long enough or does not generate sufficient force for ion channel gating. We developed a workflow using pathway similarity analysis to identify groups of residues working together to promote binding. These conformation-dependent pathways were not significantly impacted by the presence of N-linked glycans, which act primarily by interacting with the LBD bottom lobe to stabilize the closed LBD.

Fine-tuning the spike: role of the nature and topology of the glycan shield in the structure and dynamics of the SARS-CoV-2 S.

The dense glycan shield is an essential feature of the SARS-CoV-2 spike (S) architecture, key to immune evasion and to the activation of the prefusion conformation. Recent studies indicate that the occupancy and structures of the SARS-CoV-2 S glycans depend not only on the nature of the host cell, but also on the structural stability of the trimer; a point that raises important questions about the relative competence of different glycoforms. Moreover, the functional role of the glycan shield in the SARS-CoV-2 pathogenesis suggests that the evolution of the sites of glycosylation is potentially intertwined with the evolution of the protein sequence to affect optimal activity. Our results from multi-microsecond molecular dynamics simulations indicate that the type of glycosylation at N234, N165 and N343 greatly affects the stability of the receptor binding domain (RBD) open conformation, and thus its exposure and accessibility. Furthermore, our results suggest that the loss of glycosylation at N370, a newly acquired modification in the SARS-CoV-2 S glycan shield's topology, may have contributed to increase the SARS-CoV-2 infectivity as we find that N-glycosylation at N370 stabilizes the closed RBD conformation by binding a specific cleft on the RBD surface. We discuss how the absence of the N370 glycan in the SARS-CoV-2 S frees the RBD glycan binding cleft, which becomes available to bind cell-surface glycans, and potentially increases host cell surface localization.

Insights into substrate recognition and specificity for IgG by Endoglycosidase S2.

Antibodies bind foreign antigens with high affinity and specificity leading to their neutralization and/or clearance by the immune system. The conserved N-glycan on IgG has significant impact on antibody effector function, with the endoglycosidases of Streptococcus pyogenes deglycosylating the IgG to evade the immune system, a process catalyzed by the endoglycosidase EndoS2. Studies have shown that two of the four domains of EndoS2, the carbohydrate binding module (CBM) and the glycoside hydrolase (GH) domain are critical for catalytic activity. To yield structural insights into contributions of the CBM and the GH domains as well as the overall flexibility of EndoS2 to the proteins' catalytic activity, models of EndoS2-Fc complexes were generated through enhanced-sampling molecular-dynamics (MD) simulations and site-identification by ligand competitive saturation (SILCS) docking followed by reconstruction and multi-microsecond MD simulations. Modeling results predict that EndoS2 initially interacts with the IgG through its CBM followed by interactions with the GH yielding catalytically competent states. These may involve the CBM and GH of EndoS2 simultaneously interacting with either the same Fc CH2/CH3 domain or individually with the two Fc CH2/CH3 domains, with EndoS2 predicted to assume closed conformations in the former case and open conformations in the latter. Apo EndoS2 is predicted to sample both the open and closed states, suggesting that either complex can directly form following initial IgG-EndoS2 encounter. Interactions of the CBM and GH domains with the IgG are predicted to occur through both its glycan and protein regions. Simulations also predict that the Fc glycan can directly transfer from the CBM to the GH, facilitating formation of catalytically competent complexes and how the 734 to 751 loop on the CBM can facilitate extraction of the glycan away from the Fc CH2/CH3 domain. The predicted models are compared and consistent with Hydrogen/Deuterium Exchange data. In addition, the complex models are consistent with the high specificity of EndoS2 for the glycans on IgG supporting the validity of the predicted models.

Length Dependent Folding Kinetics of Alanine-Based Helical Peptides from Optimal Dimensionality Reduction.

We present a computer simulation study of helix folding in alanine homopeptides (ALA)n of length n = 5, 8, 15, and 21 residues. Based on multi-microsecond molecular dynamics simulations at room temperature, we found helix populations and relaxation times increasing from about 6% and ~2 ns for ALA5 to about 60% and ~500 ns for ALA21, and folding free energies decreasing linearly with the increasing number of residues. The helix folding was analyzed with the Optimal Dimensionality Reduction method, yielding coarse-grained kinetic models that provided a detailed representation of the folding process. The shorter peptides, ALA5 and ALA8, tended to convert directly from coil to helix, while ALA15 and ALA21 traveled through several intermediates. Coarse-grained aggregate states representing the helix, coil, and intermediates were heterogeneous, encompassing multiple peptide conformations. The folding involved multiple pathways and interesting intermediate states were present on the folding paths, with partially formed helices, turns, and compact coils. Statistically, helix initiation was favored at both termini, and the helix was most stable in the central region. Importantly, we found the presence of underlying universal local dynamics in helical peptides with correlated transitions for neighboring hydrogen bonds. Overall, the structural and dynamical parameters extracted from the trajectories are in good agreement with experimental observables, providing microscopic insights into the complex helix folding kinetics.

Computational epitope map of SARS-CoV-2 spike protein.

The primary immunological target of COVID-19 vaccines is the SARS-CoV-2 spike (S) protein. S is exposed on the viral surface and mediates viral entry into the host cell. To identify possible antibody binding sites, we performed multi-microsecond molecular dynamics simulations of a 4.1 million atom system containing a patch of viral membrane with four full-length, fully glycosylated and palmitoylated S proteins. By mapping steric accessibility, structural rigidity, sequence conservation, and generic antibody binding signatures, we recover known epitopes on S and reveal promising epitope candidates for structure-based vaccine design. We find that the extensive and inherently flexible glycan coat shields a surface area larger than expected from static structures, highlighting the importance of structural dynamics. The protective glycan shield and the high flexibility of its hinges give the stalk overall low epitope scores. Our computational epitope-mapping procedure is general and should thus prove useful for other viral envelope proteins whose structures have been characterized.

Lipid interactions of an actinoporin pore-forming oligomer

The actinoporins are cytolytic toxins produced by sea anemones. Upon encountering a membrane, preferably containing sphingomyelin, they oligomerize and insert their N-terminal helix into the membrane, forming a pore. Whether sphingomyelin is specifically recognized by the protein or simply induces phase coexistence in the membrane has been debated. Here, we perform multi-microsecond molecular dynamics simulations of an octamer of fragaceatoxin C, a member of the actinoporin family, in lipid bilayers containing either pure 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) or a 1:1 mixture of DOPC and palmitoyl sphingomyelin (PSM). The complex is highly stable in both environments, with only slight fraying of the inserted helices near their N-termini. Analyzing the structural parameters of the mixed membrane in the course of the simulation, we see signs of a phase transition for PSM in the inner leaflet of the bilayer. In both leaflets, cross-interactions between lipids of different type decrease over time. Surprisingly, the aromatic loop thought to be responsible for sphingomyelin recognition interacts more with DOPC than PSM by the end of the simulation. These results support the notion that the key membrane property that actinoporins recognize is lipid phase coexistence.

Role of D-serine Binding to the NMDA Receptor Ligand-Binding Domain

The N-methyl-D-aspartate receptor (NMDAR) is a ligand-gated ion channel that promotes synaptic plasticity critical for learning and memory. NMDAR activation is initiated by binding of two different neurotransmitter agonists to the ligand-binding domains (LBDs). Specifically, the GluN2 subunit binds glutamate, while the GluN1 subunit binds either glycine or D-serine. While the mechanisms of glycine and glutamate binding to NMDAR LBDs have already been determined, little is known about the mechanisms by which D-serine can bind to and activate the receptor. To probe binding mechanisms, we performed multi-microsecond molecular dynamics simulations of the GluN1/GluN2A LBD dimer in the presence of freely diffusing D-serine and glutamate and observed multiple binding and unbinding events. While D-serine is primarily a GluN1 agonist, we surprisingly observed that D-serine binds to both GluN1 and GluN2A NMDAR subunits. Computing interaction frequencies between LBD residues and D-serine revealed residues important for agonist binding. However, we observed that binding pathways are made up of sequences of interactions that guide D-serine into the binding site. To identify properties of these pathways, we developed a framework for analyzing binding pathways by geometric and spatial similarity, which allowed us to cluster them according to the region of the LBD with which D-serine interacts. This allowed us to identify groups of residues that work together to contribute to particular binding pathways. We determined that D-serine and glutamate bind by similar pathways and found that D-serine is able to stabilize the active LBD conformation. Lastly, we determined that N-linked glycans at physiologically relevant glycosylation sites generally do not interfere with agonist-binding pathways to modulate LBD activity. These insights highlight the versatility of D-serine as an NMDAR agonist and its unique role in receptor activation.

Thyroxine binding to type III iodothyronine deiodinase

Iodothyronine deiodinases (Dios) are important selenoproteins that control the concentration of the active thyroid hormone (TH) triiodothyronine through regioselective deiodination. The X-ray structure of a truncated monomer of Type III Dio (Dio3), which deiodinates TH inner rings through a selenocysteine (Sec) residue, revealed a thioredoxin-fold catalytic domain supplemented with an unstructured Ω-loop. Loop dynamics are driven by interactions of the conserved Trp207 with solvent in multi-microsecond molecular dynamics simulations of the Dio3 thioredoxin(Trx)-fold domain. Hydrogen bonding interactions of Glu200 with residues conserved across the Dio family anchor the loop’s N-terminus to the active site Ser-Cys-Thr-Sec sequence. A key long-lived loop conformation coincides with the opening of a cryptic pocket that accommodates thyroxine (T4) through an I⋯Se halogen bond to Sec170 and the amino acid group with a polar cleft. The Dio3-T4 complex is stabilized by an I⋯O halogen bond between an outer ring iodine and Asp211, consistent with Dio3 selectivity for inner ring deiodination. Non-conservation of residues, such as Asp211, in other Dio types in the flexible portion of the loop sequence suggests a mechanism for regioselectivity through Dio type-specific loop conformations. Cys168 is proposed to attack the selenenyl iodide intermediate to regenerate Dio3 based upon structural comparison with related Trx-fold proteins.

Markov State Model of Lassa Virus Nucleoprotein Reveals Large Structural Changes during the Trimer to Monomer Transition

Lassa virus contains a nucleoprotein (NP) that encapsulates the viral genomic RNA forming the ribonucleoprotein (RNP). The NP forms trimers that do not bind RNA, but a structure of only the NP N-terminal domain was co-crystallized with RNA bound. These structures suggested a model in which the NP forms a trimer to keep the RNA gate closed, but then is triggered to undergo a change to a form competent for RNA binding. Here, we investigate the scenario in which the trimer is disrupted to observe whether monomeric NP undergoes significant conformational changes. From multi-microsecond molecular dynamics simulations and an adaptive sampling scheme to sample the conformational space, a Markov state model (MSM) is constructed. The MSM reveals an energetically favorable conformational change, with the most significant changes occurring at the domain interface. These results support a model in which significant structural reorganization of the NP is required for RNP formation.

Solving the Gating Mechanism of the Mitochondrial β-Barrel Metabolite Channel VDAC

The Voltage-Dependent Anion Channel (VDAC) is a β-barrel protein at the mitochondrial outer membrane (MOM) serving as the major pathway for ions and respiratory metabolites between mitochondria and the cytosol. Under applied voltage, VDAC transitions from an anionic “open” state, allowing the passive diffusion of neutral and anionic metabolites, to multiple low-conducting “closed” states with reduced anion selectivity. Importantly, VDAC closed states are virtually impermeable to ATP. A potential implication of these closed states in regulating metabolite traffic across the MOM has triggered vigorous investigations on their structural nature. Still, the mechanism behind VDAC closure remains enigmatic. Here, we study VDAC1 voltage-gating by using ion-channel electrophysiology guided by ANTON2 multi-microsecond Molecular Dynamics (MD) simulations. MD simulations highlight the particular role played by conformational dynamics of a “trigger” region, which is especially impacted by voltage polarity or solution acidification. We found that a basic residue in this region is key for initiating gating transitions at neutral pH, and its mutation to the oppositely charged residue hinders channel closure and reverses selectivity. We further explore VDAC voltage-gating by combining selectivity and pH-titration experiments with Markov State Model analysis on the motions of charged residues obtained from MD simulations. We propose a VDAC gating mechanism in which all charged residues are part of an ensemble where the motion of the single basic residue at the trigger region acts as a seed coordinating distinct dynamic rearrangements of the salt-bridge complex network. Thus, VDAC closure does not imply significant conformational changes of the β-barrel. This is strikingly different from the gating mechanism of most highly ion-selective channels and thus establishes a new paradigm of β-barrel channel gating that does not involve the translocation of a voltage sensor or a major change in pore geometry.

Investigating Full-Length p53 Tetramer Dynamics with Multi-Microsecond Molecular Dynamics Simulations

Investigating Full-Length p53 Tetramer Dynamics with Multi-Microsecond Molecular Dynamics Simulations

Structures and dynamics of the novel S1/S2 protease cleavage site loop of the SARS-CoV-2 spike glycoprotein

At the end of 2019, a new highly virulent coronavirus known under the name SARS-CoV-2 emerged as a human pathogen. One key feature of SARS-CoV-2 is the presence of an enigmatic insertion in the spike glycoprotein gene representing a novel multibasic S1/S2 protease cleavage site. The proteolytic cleavage of the spike at this site is essential for viral entry into host cells. However, it has been systematically abrogated in structural studies in order to stabilize the spike in the prefusion state. In this study, multi-microsecond molecular dynamics simulations and ab initio modeling were leveraged to gain insights into the structures and dynamics of the loop containing the S1/S2 protease cleavage site. They unveiled distinct conformations, formations of short helices and interactions of the loop with neighboring glycans that could potentially regulate the accessibility of the cleavage site to proteases and its processing. In most conformations, this loop protrudes from the spike, thus representing an attractive SARS-CoV-2 specific therapeutic target.

Oxidation State Dependent Conformational Changes of HMGB1 Regulate the Formation of the CXCL12/HMGB1 Heterocomplex

High-mobility Group Box 1 (HMGB1) is an abundant protein present in all mammalian cells and involved in several processes. During inflammation or tissue damage, HMGB1 is released in the extracellular space and, depending on its redox state, can form a heterocomplex with CXCL12. The heterocomplex acts exclusively via the chemokine receptor CXCR4 enhancing leukocyte recruitment.Here, we used multi-microsecond molecular dynamics (MD) simulations to elucidate the effect of the disulfide bond on the structure and dynamics of HMGB1.The results of the MD simulations show that the presence or lack of the disulfide bond between Cys23 and Cys45 modulates the conformational space explored by HMGB1, making the reduced protein more suitable to form a complex with CXCL12.

Conformational and Dynamic Properties of Poly(ethylene oxide) in BMIM+BF4–: A Microsecond Computer Simulation Study Using ab Initio Force Fields

The behavior of polymers in complex solvents is interesting from a fundamental perspective and of practical importance from the standpoint of polymer processing. There has been recent interest in the conformational and dynamic properties of polymer in room temperature ionic liquids, with conflicting predictions from computations using models with different resolutions and conflicting results of experiments from different groups. In this work, we develop a first-principles, nonpolarizable united atom (UA) force field for a mixture of poly(ethylene oxide) (PEO) in the ionic liquid BMIM+BF4–. The UA force field is benchmarked against ab initio calculations, and the PEO atomic charges are parametrized to implicitly capture the polarization contribution to the solvation energy of a single PEO molecule in BMIM+BF4–. The UA model allows one to perform multi-microsecond molecular dynamics simulations. This is necessary because the conformational relaxation correlation times are of the order of 100 ns. The simulations predict that the radius of gyration, Rg, scales with molecular weight, Rg ∼ Mwν with ν ≈ 0.56 in the temperature range 300–600 K, consistent with experiment, seemingly in between a self-avoiding walk and an ideal chain. An examination of the snapshots of the polymer demonstrates, however, that the polymer conformations are composed of ringlike and linear segments, with ringlike parts of the chain wrapped around cations of the ionic liquid. The slow dynamics arises from the barrier to unwrapping the ringlike segments of the polymer. The mean-square displacement shows three regimes which we interpret as confinement, Zimm, and diffusive. The simulations emphasize the importance of accurate force fields and microsecond simulations in obtaining reliable results for polymers and elucidate important correlation effects for polymers in strongly interacting solvents.

Configurational Preference of the Glutamate Receptor Ligand Binding Domain Dimers

Ionotropic glutamate receptors are a family of tetrameric ion channels with functional states consisting of nonconducting, conducting, and desensitized states that are starting to become well characterized by electrophysiological and biophysical studies. However, the structure and relative energetics of these states beyond the general structure of the receptor are still not well understood. It is known that the interface between monomeric subunits of the tetramer plays a major role in distinguishing these functional states. We have used umbrella sampling and multimicrosecond molecular dynamics simulations of the GluA2 AMPA subtype glutamate receptor ligand-binding domain (LBD) dimers to characterize a natural propensity of the LBD dimers for various configurational states. Our results show a proposed desensitized conformation of the LBD dimer is a highly preferable conformation of the LBD dimer without the influence of other receptor domains or crystallographic conditions. This has been demonstrated by both free protein simulations of 5 μs duration, as well as by computed free energy difference between the active and desensitized states. At the same time, the simulations performed using the same protocols revealed that for the LBD mutant L483Y, known to lack desensitization, the postulated active state of the LBD dimer is indeed the preferred configurational state, which remained stable in the simulations. Our findings pave the path for developing more detailed hypotheses of the full receptor activation mechanism. Combined with the energetics of glutamate binding to the LBD and the energy required to open the transmembrane pore helices, our results strongly support a hypothesis that the low absolute free-energy state is the desensitized state of the intact AMPA receptor.

The Role of Unconventional Hydrogen Bonds in Determining BII Propensities in B-DNA.

An accurate understanding of DNA backbone transitions is likely to be the key for elucidating the puzzle of the intricate sequence-dependent mechanical properties that govern most of the biologically relevant functions of the double helix. One factor believed to be important in indirect recognition within protein-DNA complexes is the combined effect of two DNA backbone torsions (ε and ζ) which give rise to the well-known BI/BII conformational equilibrium. In this work we explain the sequence-dependent BII propensity observed in RpY steps (R = purine; Y = pyrimidine) at the tetranucleotide level with the help of a previously undetected C-H···O contact between atoms belonging to adjacent bases. Our results are supported by extensive multimicrosecond molecular dynamics simulations from the Ascona B-DNA Consortium, high-level quantum mechanical calculations, and data mining of the experimental structures deposited in the Protein Data Bank.

Structural Plasticity of Cholesteryl Ester Transfer Protein Assists the Lipid Transfer Activity

Cholesteryl ester transfer protein (CETP) mediates the transfer of cholesteryl esters (CEs) and triglycerides between different lipoproteins. Recent studies have shown that blocking the function of CETP can increase the level of HDL cholesterol in blood plasma and suppress the risk of cardiovascular disease. Hence, understanding the structure, dynamics, and mechanism by which CETP transfers the neutral lipids has received tremendous attention in last decade. Although the recent crystal structure has provided direct evidence of the existence of strongly bound CEs in the CETP core, very little is known about the mechanism of CE/triglyceride transfer by CETP. In this study, we explore the large scale dynamics of CETP by means of multimicrosecond molecular dynamics simulations and normal mode analysis, which provided a wealth of detailed information about the lipid transfer mechanism of CETP. Results show that the bound CEs intraconvert between bent and linear conformations in the CETP core tunnel as a consequence of the high degree of conformational flexibility of the protein. During the conformational switching, there occurred a significant reduction in hydrophobic contacts between the CEs and CETP, and a continuous tunnel traversing across the CETP long axis appeared spontaneously. Thus, our results support the recently proposed "tunnel mechanism" of CETP from cryo-EM studies for the transfer of neutral lipids between different lipoproteins. The detailed understanding obtained here could help in devising methods to prevent CETP function as a cardiovascular disease therapeutic.

Multi-Microsecond Molecular Dynamics Simulations of the HV1 Proton Channel

The voltage-gated proton channel, Hv1, is a membrane protein that plays roles in a variety of cellular processes including proton extrusion, pH regulation, and the production of reactive oxygen species. Like the voltage-gated potassium channels, it contains a voltage-sensing domain (VSD) but lacks a separate pore domain. The crystal structure of a chimera of mouse Hv1 (mHv1cc) in a non-native trimeric construct was recently reported. In order to characterize the structure and transition between the closed and open states of Hv1 in a realistic membrane environment, we used multi-microsecond molecular dynamics (MD) simulations. Specifically, we carried out MD simulations of the mouse chimera trimer and monomer in a fully hydrated POPC bilayer, as well as human Hv1 (hHv1; starting from a homology model built from mmHv1cc). We found that all three of the constructs considered (mHv1cc trimer, mHv1cc monomer, and hHv1 monomer) were stable at zero potential on the timescale of several microseconds. Additionally, the monomer hHv1 model was subjected to hyperpolarizing and depolarizing transmembrane potentials. The structure underwent only minor changes under hyperpolarizing potential. However, under a depolarizing potential, we observed a large conformational change of the protein through the outward movement of the S4 helix by approximately 8 Angstroms. This transition led to a rearrangement of the internal salt-bridge network and allowed more waters into the pore as well as an increased stability in the hydrogen bond pathways, suggesting that it is an open state. We validated our atomistic hHv1 models with a variety of experimental data, including gating charge measurements and binding of the known inhibitor, 2-GBI, at the intracellular side of the selectivity filter. Additional simulations of selected mutants under a depolarizing potential constructed from the generated open state model suggest lower proton conductivity.

Pore Formation Mechanisms of Melittin-Like Membrane-Active Peptides

Pore formation by membrane-active peptides is believed to play a key role in the mechanism of action of antimicrobial peptides. However, at present the process of spontaneous pore formation, the pore structure, and the pore lifetime in the bilayer remain unknown. This is largely due to experimental difficulties in resolving transient protein structures and dynamics in fluid-phase lipid bilayers. Here we employ a mixture of unbiased multi-microsecond molecular dynamics simulations and experimental techniques to study the process of pore formation and conductance of melittin, a pore-forming membrane active peptide from honey-bee venom, as well as synthetically evolved gain-of-function variants of melittin, which have recently been developed by William Wimley at Tulane University.Using spectroscopic techniques we show that gain-of-function variants are able to leak larger dyes from unilamellar vesicles than melittin. Like melittin gain-of-function variants result in virtually complete leakage of the dye from vesicles. However, the concentration of peptide required to achieve leakage is significantly lower for gain-of-function variants than for melittin. This suggests that these peptides form larger pores than melittin and that these pores are much more stable, remaining functional over the lifetime of the leakage experiment. Simulations of both melittin and gain-of-function variants reveal a wealth of atomic detail information about transient processes such as peptide absorption, folding, and oligomeric assembly, as well as the equilibrium structural ensemble and stability, which were verified using circular dichroism, fluorescence, and electrochemical impedance spectroscopy.

The Structure of Actin

The protein actin exists as both a monomer (G-actin) and polymer (F-actin). There have been more than 80 different crystal structures of G-actin and related homologs over the years, and at least 4 different structures/models for F-actin have been proposed using a combination of techniques. Comparison of the G-actin crystal structure with the most recent F-actin structures suggests that there are several structural rearrangements, the most significant being a relative rotation of subdomains that results in a propeller-like twist of the protein. To investigate these structures in more detail, we performed multi-microsecond molecular dynamics simulations of the muscle ADP-actin monomer (pdb 1J6Z) and the two most recent ADP F-actin structures (2ZWH and 3MFP). We find that G-actin monomer deviates quite far from its crystal structure over the course of its simulation while the F-actin protomer shows a more modest shift from its starting point. The most interesting result is that the equilibrium structures of the monomer and F-actin protomer are much more similar than previously thought and their conformational spaces show significant overlap. We will present detailed analysis of the conformational dynamics, as well as analysis of the filament conformation, twist, and structures of the barbed and pointed ends.