The Structural Basis of IKs Ion-Channel Activation: Mechanistic Insights from Molecular Simulations

Relating ion channel (iCh) structural dynamics to physiological function remains a challenge. Current experimental and computational techniques have limited ability to explore this relationship in atomistic detail over physiological timescales. A framework associating iCh structure to function is necessary for elucidating normal and disease mechanisms. We formulated a modeling schema that overcomes the limitations of current methods through applications of artificial intelligence machine learning. Using this approach, we studied molecular processes that underlie human IKs voltage-mediated gating. IKs malfunction underlies many debilitating and life-threatening diseases. Molecular components of IKs that underlie its electrophysiological function include KCNQ1 (a pore-forming tetramer) and KCNE1 (an auxiliary subunit). Simulations, using the IKs structure-function model, reproduced experimentally recorded saturation of gating-charge displacement at positive membrane voltages, two-step voltage sensor (VS) movement shown by fluorescence, iCh gating statistics, and current-voltage relationship. Mechanistic insights include the following: 1) pore energy profile determines iCh subconductance; 2) the entire protein structure, not limited to the pore, contributes to pore energy and channel subconductance; 3) interactions with KCNE1 result in two distinct VS movements, causing gating-charge saturation at positive membrane voltages and current activation delay; and 4) flexible coupling between VS and pore permits pore opening at lower VS positions, resulting in sequential gating. The new modeling approach is applicable to atomistic scale studies of other proteins on timescales of physiological function.     

A Structural Model for Apolipoprotein C-II Amyloid Fibrils: Experimental Characterization and Molecular Dynamics Simulations

The self-assembly of specific proteins to form insoluble amyloid fibrils is a characteristic feature of a number of age-related and debilitating diseases. Lipid-free human apolipoprotein C-II (apoC-II) forms characteristic amyloid fibrils and is one of several apolipoproteins that accumulate in amyloid deposits located within atherosclerotic plaques. X-ray diffraction analysis of aligned apoC-II fibrils indicated a simple cross-β-structure composed of two parallel β-sheets. Examination of apoC-II fibrils using transmission electron microscopy, scanning transmission electron microscopy, and atomic force microscopy indicated that the fibrils are flat ribbons composed of one apoC-II molecule per 4.7-Å rise of the cross-β-structure. Cross-linking results using single-cysteine substitution mutants are consistent with a parallel in-register structural model for apoC-II fibrils. Fluorescence resonance energy transfer analysis of apoC-II fibrils labeled with specific fluorophores provided distance constraints for selected donor-acceptor pairs located within the fibrils. These findings were used to develop a simple ‘letter-G-like’ β-strand-loop-β-strand model for apoC-II fibrils. Fully solvated all-atom molecular dynamics (MD) simulations showed that the model contained a stable cross-β-core with a flexible connecting loop devoid of persistent secondary structure. The time course of the MD simulations revealed that charge clusters in the fibril rearrange to minimize the effects of same-charge interactions inherent in parallel in-register models. Our structural model for apoC-II fibrils suggests that apoC-II monomers fold and self-assemble to form a stable cross-β-scaffold containing relatively unstructured connecting loops.     

Structural Basis of Type 2A von Willebrand Disease Investigated by Molecular Dynamics Simulations and Experiments

The hemostatic function of von Willebrand factor is downregulated by the metalloprotease ADAMTS13, which cleaves at a unique site normally buried in the A2 domain. Exposure of the proteolytic site is induced in the wild-type by shear stress as von Willebrand factor circulates in blood. Mutations in the A2 domain, which increase its susceptibility to cleavage, cause type 2A von Willebrand disease. In this study, molecular dynamics simulations suggest that the A2 domain unfolds under tensile force progressively through a series of steps. The simulation results also indicated that three type 2A mutations in the C-terminal half of the A2 domain, L1657I, I1628T and E1638K, destabilize the native state fold of the protein. Furthermore, all three type 2A mutations lowered in silico the tensile force necessary to undock the C-terminal helix 6 from the rest of the A2 domain, the first event in the unfolding pathway. The mutations F1520A, I1651A and A1661G were also predicted by simulations to destabilize the A2 domain and facilitate exposure of the cleavage site. Recombinant A2 domain proteins were expressed and cleavage assays were performed with the wild-type and single-point mutants. All three type 2A and two of the three predicted mutations exhibited increased rate of cleavage by ADAMTS13. These results confirm that destabilization of the helix 6 in the A2 domain facilitates exposure of the cleavage site and increases the rate of cleavage by ADAMTS13.     

Light Activation of Rhodopsin: Insights from Molecular Dynamics Simulations Guided by Solid-State NMR Distance Restraints

Structural restraints provided by solid-state NMR measurements of the metarhodopsin II intermediate are combined with molecular dynamics simulations to help visualize structural changes in the light activation of rhodopsin. Since the timescale for the formation of the metarhodopsin II intermediate (> 1 ms) is beyond that readily accessible by molecular dynamics, we use NMR distance restraints derived from ¹³C dipolar recoupling measurements to guide the simulations. The simulations yield a working model for how photoisomerization of the 11-cis retinylidene chromophore bound within the interior of rhodopsin is coupled to transmembrane helix motion and receptor activation. The mechanism of activation that emerges is that multiple switches on the extracellular (or intradiscal) side of rhodopsin trigger structural changes that converge to disrupt the ionic lock between helices H3 and H6 on the intracellular side of the receptor.     

The Possible Structural Models for Polyglutamine Aggregation: A Molecular Dynamics Simulations Study

Abstract

Huntington's disease is a neurodegenerative disorder caused by a polyglutamine (polyQ) expansion near the N-terminus of huntingtin. Previous studies have suggested that polyQ aggregation occurs only when the number of glutamine (Q) residues is more than 36-40, the disease threshold. However, the structural characteristics of polyQ nucleation in the very early stage of aggregation still remain elusive. In this study, we designed 18 simulation trials to determine the possible structural models for polyQ nucleation and aggregation with various shapes and sizes of initial β-helical structures, such as left-handed circular, right-handed rectangular, and left- and right-handed triangular. Our results show that the stability of these models significantly increases with increasing the number of rungs, while it is rather insensitive to the number of Qs in each rung. In particular, the 3-rung β-helical models are stable when they adopt the left-handed triangular and right-handed rectangular conformations due to the fact that they preserve high β-turn and β-sheet contents, respectively, during the simulation courses. Thus, we suggested that these two stable β-helical structures with at least 3 rungs might serve as the possible nucleation seeds for polyQ depending on how the structural elements of β-turn and β-sheet are sampled and preserved during the very early stage of aggregation.     

Molecular Dynamics Simulations for 1:1 Solvent Primitive Model Electrolyte Solutions

Abstract

Molecular dynamics (MD) simulations at constant temperature have been carried out for systems of 1:1 solvent primitive model (SPM) electrolyte solutions. Equilibrium thermodynamics, mean cluster size, self-diffusion coefficients, and collision frequencies were computed to examine the electrostatic effects on the structural and dynamical properties. Coherent ionic cluster motion was deduced from a cluster analysis and from the dependence of the velocity and force autocorrelation functions (FACFs). The resulting MD data for the collision frequencies and self-diffusivities of both ions and hard-spheres were shown to be in good agreement with the theoretical predictions.     

Insights into Stabilizing Forces in Amyloid Fibrils of Differing Sizes from Polarizable Molecular Dynamics Simulations

Pathological aggregation of amyloid-forming proteins is a hallmark of a number of human diseases, including Alzheimer's, type 2 diabetes, Parkinson's, and more. Despite having very different primary amino acid sequences, these amyloid proteins form similar supramolecular, fibril structures that are highly resilient to physical and chemical denaturation. To better understand the structural stability of disease-related amyloids and to gain a greater understanding of factors that stabilize functional amyloid assemblies, insights into tertiary and quaternary interactions are needed. We performed molecular dynamics simulations on human tau, amyloid-β, and islet amyloid polypeptide fibrils to determine key physicochemical properties that give rise to their unique characteristics and fibril structures. These simulations are the first of their kind in employing a polarizable force field to explore properties of local electric fields on dipole properties and other electrostatic forces that contribute to amyloid stability. Across these different amyloid fibrils, we focused on how the underlying forces stabilize fibrils to elucidate the driving forces behind the protein aggregation. The polarizable model allows for an investigation of how side-chain dipole moments, properties of structured water molecules in the fibril core, and the local environment around salt bridges contribute to the formation of interfaces essential for fibril stability. By systematically studying three amyloidogenic proteins of various fibril sizes for key structural properties and stabilizing forces, we shed light on properties of amyloid structures related to both diseased and functional states at the atomistic level.     

A Kinetic Model of Trp-Cage Folding from Multiple Biased Molecular Dynamics Simulations

Trp-cage is a designed 20-residue polypeptide that, in spite of its size, shares several features with larger globular proteins. Although the system has been intensively investigated experimentally and theoretically, its folding mechanism is not yet fully understood. Indeed, some experiments suggest a two-state behavior, while others point to the presence of intermediates. In this work we show that the results of a bias-exchange metadynamics simulation can be used for constructing a detailed thermodynamic and kinetic model of the system. The model, although constructed from a biased simulation, has a quality similar to those extracted from the analysis of long unbiased molecular dynamics trajectories. This is demonstrated by a careful benchmark of the approach on a smaller system, the solvated Ace-Ala₃-Nme peptide. For the Trp-cage folding, the model predicts that the relaxation time of 3100 ns observed experimentally is due to the presence of a compact molten globule-like conformation. This state has an occupancy of only 3% at 300 K, but acts as a kinetic trap. Instead, non-compact structures relax to the folded state on the sub-microsecond timescale. The model also predicts the presence of a state at of 4.4 Å from the NMR structure in which the Trp strongly interacts with Pro12. This state can explain the abnormal temperature dependence of the and chemical shifts. The structures of the two most stable misfolded intermediates are in agreement with NMR experiments on the unfolded protein. Our work shows that, using biased molecular dynamics trajectories, it is possible to construct a model describing in detail the Trp-cage folding kinetics and thermodynamics in agreement with experimental data.     

Inherent Dynamics of Head Domain Correlates with ATP-Recognition of P2X4 Receptors: Insights Gained from Molecular Simulations

P2X receptors are ATP-gated ion channels involved in many physiological functions, and determination of ATP-recognition (AR) of P2X receptors will promote the development of new therapeutic agents for pain, inflammation, bladder dysfunction and osteoporosis. Recent crystal structures of the zebrafish P2X4 (zfP2X4) receptor reveal a large ATP-binding pocket (ABP) located at the subunit interface of zfP2X4 receptors, which is occupied by a conspicuous cluster of basic residues to recognize triphosphate moiety of ATP. Using the engineered affinity labeling and molecular modeling, at least three sites (S1, S2 and S3) within ABP have been identified that are able to recognize the adenine ring of ATP, implying the existence of at least three distinct AR modes in ABP. The open crystal structure of zfP2X4 confirms one of three AR modes (named AR1), in which the adenine ring of ATP is buried into site S1 while the triphosphate moiety interacts with clustered basic residues. Why architecture of ABP favors AR1 not the other two AR modes still remains unexplored. Here, we examine the potential role of inherent dynamics of head domain, a domain involved in ABP formation, in AR determinant of P2X4 receptors. In silico docking and binding free energy calculation revealed comparable characters of three distinct AR modes. Inherent dynamics of head domain, especially the downward motion favors the preference of ABP for AR1 rather than AR2 and AR3. Along with the downward motion of head domain, the closing movement of loop_139–146 and loop_169–183, and structural rearrangements of K70, K72, R298 and R143 enabled ABP to discriminate AR1 from other AR modes. Our observations suggest the essential role of head domain dynamics in determining AR of P2X4 receptors, allowing evaluation of new strategies aimed at developing specific blockers/allosteric modulators by preventing the dynamics of head domain associated with both AR and channel activation of P2X4 receptors.     

Correlation Spectroscopy and Molecular Dynamics Simulations to Study the Structural Features of Proteins

In this work, we used a combination of fluorescence correlation spectroscopy (FCS) and molecular dynamics (MD) simulation methodologies to acquire structural information on pH-induced unfolding of the maltotriose-binding protein from Thermus thermophilus (MalE2). FCS has emerged as a powerful technique for characterizing the dynamics of molecules and it is, in fact, used to study molecular diffusion on timescale of microsecond and longer. Our results showed that keeping temperature constant, the protein diffusion coefficient decreased from 84±4 µm²/s to 44±3 µm²/s when pH was changed from 7.0 to 4.0. An even more marked decrease of the MalE2 diffusion coefficient (31±3 µm²/s) was registered when pH was raised from 7.0 to 10.0. According to the size of MalE2 (a monomeric protein with a molecular weight of 43 kDa) as well as of its globular native shape, the values of 44 µm²/s and 31 µm²/s could be ascribed to deformations of the protein structure, which enhances its propensity to form aggregates at extreme pH values. The obtained fluorescence correlation data, corroborated by circular dichroism, fluorescence emission and light-scattering experiments, are discussed together with the MD simulations results.     

Insights into the N-Sulfation Mechanism: Molecular Dynamics Simulations of the N-Sulfotransferase Domain of Ndst1 and Mutants

Sulfation patterns along glycosaminoglycan (GAG) chains dictate their functional role. The N-deacetylase N-sulfotransferase family (NDST) catalyzes the initial downstream modification of heparan sulfate and heparin chains by removing acetyl groups from subsets of N-acetylglucosamine units and, subsequently, sulfating the residual free amino groups. These enzymes transfer the sulfuryl group from 3′-phosphoadenosine-5′-phosphosulfate (PAPS), yielding sulfated sugar chains and 3′-phosphoadenosine-5′-phosphate (PAP). For the N-sulfotransferase domain of NDST1, Lys833 has been implicated to play a role in holding the substrate glycan moiety close to the PAPS cofactor. Additionally, Lys833 together with His716 interact with the sulfonate group, stabilizing the transition state. Such a role seems to be shared by Lys614 through donation of a proton to the bridging oxygen of the cofactor, thereby acting as a catalytic acid. However, the relevance of these boundary residues at the hydrophobic cleft is still unclear. Moreover, whether Lys833, His716 and Lys614 play a role in both glycan recognition and glycan sulfation remains elusive. In this study we evaluate the contribution of NDST mutants (Lys833, His716 and Lys614) to dynamical effects during sulfate transfer using comprehensive combined docking and essential dynamics. In addition, the binding location of the glycan moiety, PAPS and PAP within the active site of NDST1 throughout the sulfate transfer were determined by intermediate state analysis. Furthermore, NDST1 mutants unveiled Lys833 as vital for both the glycan binding and subsequent N-sulfotransferase activity of NDST1.     

Structural and Dynamical Properties of a Full-length HIV-1 Integrase: Molecular Dynamics Simulations

Abstract

The structural and dynamical properties of the complete full-length structure of HIV-1 integrase were investigated using Molecular Dynamics approach. Simulations were carried out for the three systems, core domain only (CORE), full-length structure without (FULL) and with a Mg²⁺ (FULL+ION) in its active site, aimed to investigate the difference in the molecular properties of the full-length models due to their different construction procedures as well as the effects of the two ends, C- and N-terminal, on those properties in the core domain. The full-length structure was prepared from the two experimental structures of two-domain fragment. The following properties were observed to differ significantly from the previous reports: (i) relative topology formed by an angle between the three domains; (ii) the cavity size defined by the catalytic triad, Asp64, Asp116, and Glul52; (iii) distances and solvation of the Mg²⁺; and (iv) conformation of the catalytic residues. In addition, the presence of the two terminal domains decreases the mobility of the central core domain significantly.     

Molecular Dynamics Simulations to Provide Insights into Epitopes Coupled to the Soluble and Membrane-Bound MHC-II Complexes

Epitope recognition by major histocompatibility complex II (MHC-II) is essential for the activation of immunological responses to infectious diseases. Several studies have demonstrated that this molecular event takes place in the MHC-II peptide-binding groove constituted by the α and β light chains of the heterodimer. This MHC-II peptide-binding groove has several pockets (P1-P11) involved in peptide recognition and complex stabilization that have been probed through crystallographic experiments and in silico calculations. However, most of these theoretical calculations have been performed without taking into consideration the heavy chains, which could generate misleading information about conformational mobility both in water and in the membrane environment. Therefore, in absence of structural information about the difference in the conformational changes between the peptide-free and peptide-bound states (pMHC-II) when the system is soluble in an aqueous environment or non-covalently bound to a cell membrane, as the physiological environment for MHC-II is. In this study, we explored the mechanistic basis of these MHC-II components using molecular dynamics (MD) simulations in which MHC-II was previously co-crystallized with a small epitope (P₇) or coupled by docking procedures to a large (P₂₂) epitope. These MD simulations were performed at 310 K over 100 ns for the water-soluble (MHC-II_w, MHC-II-P_7w, and MHC-II-P_22w) and 150 ns for the membrane-bound species (MHC-II_m, MHC-II-P_7m, and MHC-II-P_22m). Our results reveal that despite the different epitope sizes and MD simulation environments, both peptides are stabilized primarily by residues lining P1, P4, and P6-7, and similar noncovalent intermolecular energies were observed for the soluble and membrane-bound complexes. However, there were remarkably differences in the conformational mobility and intramolecular energies upon complex formation, causing some differences with respect to how the two peptides are stabilized in the peptide-binding groove.     

The Fip35 WW Domain Folds with Structural and Mechanistic Heterogeneity in Molecular Dynamics Simulations

We describe molecular dynamics simulations resulting in the folding the Fip35 Hpin1 WW domain. The simulations were run on a distributed set of graphics processors, which are capable of providing up to two orders of magnitude faster computation than conventional processors. Using the Folding@home distributed computing system, we generated thousands of independent trajectories in an implicit solvent model, totaling over 2.73 ms of simulations. A small number of these trajectories folded; the folding proceeded along several distinct routes and the system folded into two distinct three-stranded β-sheet conformations, showing that the folding mechanism of this system is distinctly heterogeneous.     

The Role of Molecular Dynamics Simulations for the Study of Slow Dynamics

Abstract

A two step strategy is proposed to study dynamical properties of a physical system much slower than the time scales accessible by molecular dynamics simulations. The strategy is applied to investigate the slow dynamics of supercooled liquids.     

Long-Distance Communication in the HDV Ribozyme: Insights from Molecular Dynamics and Experiments

The hepatitis delta virus ribozyme is a small, self-cleaving RNA with a compact tertiary structure and buried active site that is important in the life cycle of the virus. The ribozyme's function in nature is to cleave an internal phosphodiester bond and linearize concatemers during rolling circle replication. Crystal structures of the ribozyme have been solved in both pre-cleaved and post-cleaved (product) forms and reveal an intricate network of interactions that conspire to catalyze bond cleavage. In addition, extensive biochemical studies have been performed to work out a mechanism for bond cleavage in which C75 and a magnesium ion catalyze the reaction by general acid-base chemistry. One issue that has remained unclear in this ribozyme and in other ribozymes is the nature of long-distance communication between peripheral regions of the RNA and the buried active site. We performed molecular dynamics simulations on the hepatitis delta virus ribozyme in the product form and assessed communication between a distal structural portion of the ribozyme—the protonated C41 base triple—and the active site containing the critical C75. We varied the ionization state of C41 in both the wild type and a C41 double mutant variant and determined the impact on the active site. In all four cases, effects at the active site observed in the simulations agree with experimental studies on ribozyme activity. Overall, these studies indicate that small functional RNAs have the potential to communicate interactions over long distances and that wild-type RNAs may have evolved ways to prevent such interactions from interfering with catalysis.     

Effects of the Pathogenic Mutation A117V and the Protective Mutation H111S on the Folding and Aggregation of PrP106-126: Insights from Replica Exchange Molecular Dynamics Simulations

The fragment 106-126 of prion protein exhibits similar properties to full-length prion. Experiments have shown that the A117V mutation enhances the aggregation of PrP106-126, while the H111S mutation abolishes the assembly. However, the mechanism of the change in the aggregation behavior of PrP106-126 upon the two mutations is not fully understood. In this study, replica exchange molecular dynamics simulations were performed to investigate the conformational ensemble of the WT PrP106-126 and its two mutants A117V and H111S. The obtained results indicate that the three species are all intrinsically disordered but they have distinct morphological differences. The A117V mutant has a higher propensity to form β-hairpin structures than the WT, while the H111S mutant has a higher population of helical structures. Furthermore, the A117V mutation increases the hydrophobic solvent accessible surface areas of PrP106-126 and the H111S mutation reduces the exposure of hydrophobic residues. It can be concluded that the difference in populations of β-hairpin structures and the change of hydrophobic solvent accessible areas may induce the different aggregation behaviors of the A117V and the H111S mutated PrP106-126. Understanding why the two mutations have contrary effects on the aggregation of PrP106-126 is very meaningful for further elucidation of the mechanism underlying aggregation and design of inhibitor against aggregation process.     

Full Length Vpu from HIV-1: Combining Molecular Dynamics Simulations with NMR Spectroscopy

Abstract

Based on structures made available by solution NMR, molecular models of the protein Vpu from HIV-1 were built and refined by 6 ns MD simulations in a fully hydrated lipid bilayer. Vpu is an 81 amino acid type I integral membrane protein encoded by the human immunodeficiency virus type-1 (HIV-1) and closely related simian immunodeficiency viruses (SIVs). Its role is to amplify viral release. Upon phosphorylation, the cytoplasmic domain adopts a more compact shape with helices 2 and 3 becoming almost parallel to each other. A loss of helicity for several residues belonging to the helices adjacent to both ends of the loop region containing serines 53 and 57 is observed. A fourth helix, present in one of the NMR-based structures of the cytoplasmic domain and located near the C-terminus, is lost upon phosphorylation.