Explicit Solvent Molecular Dynamics Research Articles

Free energy predictions of ligand binding to an α-helix using steered molecular dynamics and umbrella sampling simulations.

Free energy prediction of ligand binding to macromolecules using explicit solvent molecular dynamics (MD) simulations is computationally very expensive. Recently, we reported a linear correlation between the binding free energy obtained via umbrella sampling (US) versus the rupture force from steered molecular dynamics (SMD) simulations for epigallocatechin-3-gallate (EGCG) binding to α-helical-rich keratin. This linear correlation suggests a potential route for fast free energy predictions using SMD alone. In this work, the generality of the linear correlation is further tested for several ligands interacting with the α-helical motif of keratin. These molecules have significantly varying properties, i.e., octanol/water partition coefficient (log P), and/or overall charges (oleic acid, catechin, Fe(2+), citric acid, hydrogen citrate, dihydrogen citrate, and citrate). Using the constant loading rate of our previous study of the keratin-EGCG system, we observe that the linear correlation for keratin-EGCG can be extended to other uncharged molecules where interactions are governed by hydrogen bonds and/or a combination of hydrogen bonds and hydrophobic forces. For molecules where interactions with the keratin helix are governed primarily by electrostatics between charged molecules, a second, alternative linear correlation model is derived. While further investigations are needed to expand the molecular space and build a fully predictive model, the current approach represents a promising methodology for fast free energy predictions based on short SMD simulations (requiring picoseconds to nanoseconds of sampling) for defined biomolecular systems.

Mechanism of autophosphorylation of mycobacterial PknB explored by molecular dynamics simulations.

Mycobacterial Ser/Thr kinase, PknB, is essential for the growth of the pathogen. Unphosphorylated PknB is catalytically inactive, and its activation requires autophosphorylation of Thr residues on the activation loop. Autophosphorylation can in principle take place via two distinct mechanisms. Intermolecular trans autophosphorylation involves dimerization and phosphorylation of the activation loop of one chain in the catalytic pocket of the other chain. On the other hand, intramolecular cis autophosphorylation involves phosphorylation of the activation loop of the kinases in its own catalytic pocket within a monomer. On the basis of the crystal structure of PknB in the front-to-front dimeric form, it is currently believed that activation of PknB involves trans autophosphorylation. However, because of the lack of coordinates of the activation loop in the crystal structures, atomic details of the conformational changes associated with activation are yet to be deciphered. Therefore, to understand the conformational transitions associated with activation via autophosphorylation, a series of explicit solvent molecular dynamics simulations with a duration of 1 μs have been performed on each of the phosphorylated and nonphosphorylated forms of the PknB catalytic domain in monomeric and dimeric states. Simulations on phosphorylated PknB revealed a differential network of crucial electrostatic and hydrophobic residues that stabilize the phosphorylated form in the active conformation. Interestingly, in our simulations on nonphosphorylated monomers, the activation loop was observed to fold into its own active site, thereby opening the novel possibility of activation through intramolecular cis autophosphorylation. Thus, our simulations suggest that autophosphorylation of PknB might also involve cis initiation followed by trans amplification as reported for other eukaryotic kinases based on recent reaction kinetics studies.

Possible dynamically gated conductance along heme wires in bacterial multiheme cytochromes.

The staggered cross decaheme configuration of electron transfer cofactors in the outer-membrane cytochrome MtrF serves as a prototype for conformationally gated multiheme electron transport. Derived from the bacterium Shewanella oneidensis, the staggered cross configuration reveals intersecting c-type octaheme and tetraheme "wires" containing thermodynamic "hills" and "valleys" (Proc. Natl. Acad. Sci. U. S. A. 2014, 11, 611-616), suggesting that the protein structure may include a dynamical mechanism for conductance and pathway switching depending on enzymatic functional need. Here, we applied classical molecular and statistical mechanics calculations of large-amplitude protein dynamics in MtrF, to address its potential to modulate pathway conductance, including assessment of the effect of the total charge state. Explicit solvent molecular dynamics simulations of fully oxidized and fully reduced MtrF showed that the slowest mode of collective decaheme motion is 90% similar between the oxidized and reduced states and consists primarily of interheme separation with minor rotational contributions. The frequency of this motion is 1.7 × 10(7) s(-1), both for fully oxidized and fully reduced MtrF, slower than the downhill electron transfer rates between stacked heme pairs at the octaheme termini and faster than the electron transfer rates between parallel hemes in the tetraheme chain. This implies that MtrF uses slow conformational fluctuations to modulate electron flow along the octaheme pathway, apparently for the purpose of increasing the residence time of electrons on lowest potential hemes 4 and 9. This apparent gating mechanism should increase the success rate of electron transfer from MtrF to low potential environmental acceptors via these two solvent-exposed hemes.

Solvent friction effects propagate over the entire protein molecule through low-frequency collective modes.

Protein solvation dynamics has been investigated using atom-dependent Langevin friction coefficients derived directly from molecular dynamics (MD) simulations. To determine the effect of solvation on the atomic friction coefficients, solution and vacuum MD simulations were performed for lysozyme and staphylococcal nuclease and analyzed by Langevin mode analysis. The coefficients thus derived are roughly correlated with the atomic solvent-accessible surface area (ASA), as expected from the fact that friction occurs as the result of collisions with solvent molecules. However, a considerable number of atoms with higher friction coefficients are found inside the core region. Hence, the influence of solvent friction propagates into the protein core. The internal coefficients have large contributions from the low-frequency modes, yielding a simple picture of the surface-to-core long-range damping via solvation governed by collective low-frequency modes. To make use of these findings in implicit-solvent modeling, we compare the all-atom friction results with those obtained using Langevin dynamics (LD) with two empirical representations: the constant-friction and the ASA-dependent (Pastor-Karplus) friction models. The constant-friction model overestimates the core and underestimates the surface damping whereas the ASA-dependent friction model, which damps protein atoms only on the solvent-accessible surface, reproduces well the friction coefficients for both the surface and core regions observed in the explicit-solvent MD simulations. Therefore, in LD simulation, the solvent friction coefficients should be imposed only on the protein surface.

Validating Solution Ensembles from Molecular Dynamics Simulation by Wide-Angle X-ray Scattering Data

Wide-angle x-ray scattering (WAXS) experiments of biomolecules in solution have become increasingly popular because of technical advances in light sources and detectors. However, the structural interpretation of WAXS profiles is problematic, partly because accurate calculations of WAXS profiles from structural models have remained challenging. In this work, we present the calculation of WAXS profiles from explicit-solvent molecular dynamics (MD) simulations of five different proteins. Using only a single fitting parameter that accounts for experimental uncertainties because of the buffer subtraction and dark currents, we find excellent agreement to experimental profiles both at small and wide angles. Because explicit solvation eliminates free parameters associated with the solvation layer or the excluded solvent, which would require fitting to experimental data, we minimize the risk of overfitting. We further find that the influence from water models and protein force fields on calculated profiles are insignificant up to q≈15nm−1. Using a series of simulations that allow increasing flexibility of the proteins, we show that incorporating thermal fluctuations into the calculations significantly improves agreement with experimental data, demonstrating the importance of protein dynamics in the interpretation of WAXS profiles. In addition, free MD simulations up to one microsecond suggest that the calculated profiles are highly sensitive with respect to minor conformational rearrangements of proteins, such as an increased flexibility of a loop or an increase of the radius of gyration by < 1%. The present study suggests that quantitative comparison between MD simulations and experimental WAXS profiles emerges as an accurate tool to validate solution ensembles of biomolecules.

Enrichment of chemical libraries docked to protein conformational ensembles and application to aldehyde dehydrogenase 2.

Molecular recognition is a complex process that involves a large ensemble of structures of the receptor and ligand. Yet, most structure-based virtual screening is carried out on a single structure typically from X-ray crystallography. Explicit-solvent molecular dynamics (MD) simulations offer an opportunity to sample multiple conformational states of a protein. Here we evaluate our recently developed scoring method SVMSP in its ability to enrich chemical libraries docked to MD structures of seven proteins from the Directory of Useful Decoys (DUD). SVMSP is a target-specific rescoring method that combines machine learning with statistical potentials. We find that enrichment power as measured by the area under the ROC curve (ROC-AUC) is not affected by increasing the number of MD structures. Among individual MD snapshots, many exhibited enrichment that was significantly better than the crystal structure, but no correlation between enrichment and structural deviation from crystal structure was found. We followed an innovative approach by training SVMSP scoring models using MD structures (SVMSPMD). The resulting models were applied to two difficult cases (p38 and CDK2) for which enrichment was not better than random. We found remarkable increase in enrichment power, particularly for p38, where the ROC-AUC increased by 0.30 to 0.85. Finally, we explored approaches for a priori identification of MD snapshots with high enrichment power from an MD simulation in the absence of active compounds. We found that the use of randomly selected compounds docked to the target of interest using SVMSP led to notable enrichment for EGFR and Src MD snapshots. SVMSP rescoring of protein–compound MD structures was applied for the search of small-molecule inhibitors of the mitochondrial enzyme aldehyde dehydrogenase 2 (ALDH2). Rank-ordering of a commercial library of 50 000 compounds docked to MD structures of ALDH2 led to five small-molecule inhibitors. Four compounds had IC50s below 5 μM. These compounds serve as leads for the design and synthesis of more potent and selective ALDH2 inhibitors.

Influence of Sequence and Covalent Modifications on Yeast tRNA Dynamics.

Modified nucleotides are prevalent in tRNA. Experimental studies reveal that these covalent modifications play an important role in tuning tRNA function. In this study, molecular dynamics (MD) simulations were used to investigate how modifications alter tRNA dynamics. The X-ray crystal structures of tRNA(Asp), tRNA(Phe), and tRNA(iMet), both with and without modifications, were used as initial structures for 333 ns explicit solvent MD simulations with AMBER. For each tRNA molecule, three independent trajectory calculations were performed, giving an aggregate of 6 μs of total MD across six molecules. The global root-mean-square deviations (RMSD) of atomic positions show that modifications only introduce significant rigidity to the global structure of tRNA(Phe). Interestingly, RMSDs of the anticodon stem-loop (ASL) suggest that modified tRNA has a more rigid structure compared to the unmodified tRNA in this domain. The anticodon RMSDs of the modified tRNAs, however, are higher than those of corresponding unmodified tRNAs. These findings suggest that the rigidity of the anticodon stem-loop is finely tuned by modifications, where rigidity in the anticodon arm is essential for tRNA translocation in the ribosome, and flexibility of the anticodon is important for codon recognition. Sugar pucker and water residence time of pseudouridines in modified tRNAs and corresponding uridines in unmodified tRNAs were assessed, and the results reinforce that pseudouridine favors the 3′-endo conformation and has a higher tendency to interact with water. Principal component analysis (PCA) was used to examine correlated motions in tRNA. Additionally, covariance overlaps of PCAs were compared for trajectories of the same molecule and between trajectories of modified and unmodified tRNAs. The comparison suggests that modifications alter the correlated motions. For the anticodon bases, the extent of stacking was compared between modified and unmodified molecules, and only unmodified tRNA(Asp) has significantly higher percentage of stacking time. Overall, the simulations reveal that the effect of covalent modification on tRNA dynamics is not simple, with modifications increasing flexibility in some regions of the structure and increasing rigidity in other regions.

Structural and Energetic Insight into the Cross-Seeding Amyloid Assemblies of Human IAPP and Rat IAPP

The misfolding and aggregation of human islet amyloid polypeptide (hIAPP or amylin) into small oligomers and large amyloid fibrils is believed to be responsible for the dysfunction and death of pancreatic β-cells in diabetes type II. However, rat IAPP (rIAPP), which differs from the hIAPP by only 6 of 37 residues, lacks the ability to form amyloid fibrils and to induce cell death. Little is known about the cross-sequence interactions and cross-seeding structures between hIAPP and rIAPP peptides. Herein using explicit-solvent molecular dynamics (MD) simulations, we modeled and simulated different heteroassemblies formed by the amyloidogenic hIAPP and the nonamyloidogenic rIAPP peptides. Simulations showed that the U-shaped hIAPP monomer and oligomers can interact with conformationally similar rIAPP to form stable complexes and to coassemble into heterogeneous structures. Stable heterointeractions between hIAPP and rIAPP were shown to arise from hydrophobic contacts and hydrogen bonds at the interface, particularly at N- and C-terminal β-sheet regions. Because of the enhanced interpeptide interactions at the interface, upon binding to hIAPP oligomers, the β-sheet population of rIAPP was greatly increased as compared to that of rIAPP alone. More importantly, the conformational energies of rIAPP monomers at the bound state were observed to be consistently higher than those of rIAPP monomers at the unbound state. However, rIAPP monomers enable one to adopt different conformations and follow different pathways for associating with hIAPP from the high energy of the bound state to the low energy of the unbound state, without encountering any large and abrupt energy barrier. In parallel, AFM study of cross-aggregation of hIAPP and rIAPP provided additional evidence that hIAPP can seed with rIAPP to form hybrid fibrils at all concentrations similar to pure hIAPP fibrils. This work demonstrates the existence of cross-interactions between the two different IAPP peptides, which provides an improved fundamental understanding of the cross-seeding of different amyloid sequences toward amyloid aggregation and toxicity mechanisms.

Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange.

We present a new computational approach for constant pH simulations in explicit solvent based on the combination of the enveloping distribution sampling (EDS) and Hamiltonian replica exchange (HREX) methods. Unlike constant pH methods based on variable and continuous charge models, our method is based on discrete protonation states. EDS generates a hybrid Hamiltonian of different protonation states. A smoothness parameter s is used to control the heights of energy barriers of the hybrid-state energy landscape. A small s value facilitates state transitions by lowering energy barriers. Replica exchange between EDS potentials with different s values allows us to readily obtain a thermodynamically accurate ensemble of multiple protonation states with frequent state transitions. The analysis is performed with an ensemble obtained from an EDS Hamiltonian without smoothing, s = ∞, which strictly follows the minimum energy surface of the end states. The accuracy and efficiency of this method is tested on aspartic acid, lysine, and glutamic acid, which have two protonation states, a histidine with three states, a four-residue peptide with four states, and snake cardiotoxin with eight states. The pKa values estimated with the EDS-HREX method agree well with the experimental pKa values. The mean absolute errors of small benchmark systems range from 0.03 to 0.17 pKa units, and those of three titratable groups of snake cardiotoxin range from 0.2 to 1.6 pKa units. This study demonstrates that EDS-HREX is a potent theoretical framework, which gives the correct description of multiple protonation states and good calculated pKa values.

Interactions of a water-soluble fullerene derivative with amyloid-β protofibrils: dynamics, binding mechanism, and the resulting salt-bridge disruption.

Alzheimer's disease (AD) is associated with the pathological self-assembly of amyloid-β (Aβ) peptides into β-sheet-rich oligomers and insoluble amyloid fibrils. Experimental studies reported that 1,2-(dimethoxymethano)fullerene (DMF), a water-soluble fullerene derivative, inhibits strongly Aβ peptide aggregation at the early stage. However, the interaction and binding mechanisms are not well understood. In this study, we have investigated the detailed interaction of a DMF molecule with a fibrillar hexamer of full-length Aβ42 and the resulting structural alterations by performing multiple all-atom explicit solvent molecular dynamics (MD) simulations. Starting from different initial states with a minimum distance of 2 nm between the DMF and the Aβ protofibril, our MD simulations show that the DMF binds to the Aβ protofibril via both slow and fast binding processes. Three dominant binding sites are identified: the central hydrophobic core (CHC) site (17LVFFA21), the turn site (27NKGAI31), and the C-terminal β-sheet site consisting of the smallest side-chain residue glycine and hydrophobic residues (31IIGLMVGGVVI41). Binding energy analyses reveal the importance of π-stacking interactions, especially in the CHC site, hydrophobic interactions, and curvature matching. Strikingly, we find that the binding of DMF to the turn region can disrupt the D23-K28 salt-bridge that is important for the Aβ fibrillation. These results provide molecular insight into the binding mechanism of fullerene to Aβ protofibrils and offer new routes for the therapeutic drug design using fullerene derivatives against AD.

Two structural scenarios for protein stabilization by PEG.

PEGylation, or addition of poly(ethylene glycol) chains to proteins, is widely used to improve delivery in pharmaceutical applications. Recent studies suggest that stabilization of a protein by PEG, and hence its proteolytic degradability, is sequence-dependent and requires only short PEG chains. Here we connect stabilization by short PEG chains directly to the structural dynamics of the protein and PEG chain. We measured the stability of human Pin1 WW domain with PEG-4 at asparagine 19 for a full mutant cycle at two positions thought to influence PEG-protein interaction: Ser16Ala and Tyr23Phe. We then performed explicit solvent molecular dynamics simulations on all PEGylated and PEG-free mutants. The mutant cycle yields a nonadditive stabilization effect where the pseudo-wild type and double mutant are more stabilized relative to unPEGylated proteins than are the two single mutants. The simulation reveals why: the double mutant suffers loss of β-sheet structure, which PEGylation restores even though the PEG extends as a coil into the solvent. In contrast, in one of the single mutants, PEG preferentially interacts with the protein surface while disrupting the interactions of its asparagine host with a nearby methionine side chain. Thus, PEG attachment can stabilize a protein differentially depending on the local sequence, and either by interacting with the surface or by extending into the solvent. A simulation with PEG-45 attached to asparagine 19 shows that PEG even can do both in the same context.

Are current atomistic force fields accurate enough to study proteins in crowded environments?

The high concentration of macromolecules in the crowded cellular interior influences different thermodynamic and kinetic properties of proteins, including their structural stabilities, intermolecular binding affinities and enzymatic rates. Moreover, various structural biology methods, such as NMR or different spectroscopies, typically involve samples with relatively high protein concentration. Due to large sampling requirements, however, the accuracy of classical molecular dynamics (MD) simulations in capturing protein behavior at high concentration still remains largely untested. Here, we use explicit-solvent MD simulations and a total of 6.4 µs of simulated time to study wild-type (folded) and oxidatively damaged (unfolded) forms of villin headpiece at 6 mM and 9.2 mM protein concentration. We first perform an exhaustive set of simulations with multiple protein molecules in the simulation box using GROMOS 45a3 and 54a7 force fields together with different types of electrostatics treatment and solution ionic strengths. Surprisingly, the two villin headpiece variants exhibit similar aggregation behavior, despite the fact that their estimated aggregation propensities markedly differ. Importantly, regardless of the simulation protocol applied, wild-type villin headpiece consistently aggregates even under conditions at which it is experimentally known to be soluble. We demonstrate that aggregation is accompanied by a large decrease in the total potential energy, with not only hydrophobic, but also polar residues and backbone contributing substantially. The same effect is directly observed for two other major atomistic force fields (AMBER99SB-ILDN and CHARMM22-CMAP) as well as indirectly shown for additional two (AMBER94, OPLS-AAL), and is possibly due to a general overestimation of the potential energy of protein-protein interactions at the expense of water-water and water-protein interactions. Overall, our results suggest that current MD force fields may distort the picture of protein behavior in biologically relevant crowded environments.

Molecular Dynamics Simulations of Nucleic Acids. From Tetranucleotides to the Ribosome.

We present a brief overview of explicit solvent molecular dynamics (MD) simulations of nucleic acids. We explain physical chemistry limitations of the simulations, namely, the molecular mechanics (MM) force field (FF) approximation and limited time scale. Further, we discuss relations and differences between simulations and experiments, compare standard and enhanced sampling simulations, discuss the role of starting structures, comment on different versions of nucleic acid FFs, and relate MM computations with contemporary quantum chemistry. Despite its limitations, we show that MD is a powerful technique for studying the structural dynamics of nucleic acids with a fast growing potential that substantially complements experimental results and aids their interpretation.

QM Computations on Complete Nucleic Acids Building Blocks: Analysis of the Sarcin-Ricin RNA Motif Using DFT-D3, HF-3c, PM6-D3H, and MM Approaches.

A set of conformations obtained from explicit solvent molecular dynamics (MD) simulations of the Sarcin-Ricin internal loop (SRL) RNA motif is investigated using quantum mechanical (QM, TPSS-D3/def2-TZVP DFT-D3) and molecular mechanics (MM, AMBER parm99bsc0+χol3 force field) methods. Solvent effects are approximated using implicit solvent methods (COSMO for DFT-D3; GB and PB for MM). Large-scale DFT-D3 optimizations of the full 11-nucleotide motif are compared to MM results and reveal a higher flexibility of DFT-D3 over the MM in the optimization procedure. Conformational energies of the SRL motif expose significant differences in the DFT-D3 and MM energy descriptions that explain difficulties in MD simulations of the SRL motif. The TPSS-D3 data are in excellent agreement with results obtained by the hybrid functionals PW6B95-D3 and M06-2X. Computationally more efficient methods such as PM6-D3H and HF-3c show promising but partly inconsistent results. It is demonstrated that large-scale DFT-D3 computations on complete nucleic acids building blocks are a viable tool to complement the picture obtained from MD simulations and can be used as benchmarks for faster computational methods. Methodological challenges of large-scale QM computations on nucleic acids such as missing solvent-solute interactions and the truncation of the studied systems are discussed.

Mesoscale properties of clay aggregates from potential of mean force representation of interactions between nanoplatelets

Face-to-face and edge-to-edge free energy interactions of Wyoming Na-montmorillonite platelets were studied by calculating potential of mean force along their center to center reaction coordinate using explicit solvent (i.e., water) molecular dynamics and free energy perturbation methods. Using a series of configurations, the Gay-Berne potential was parametrized and used to examine the meso-scale aggregation and properties of platelets that are initially random oriented under isothermal-isobaric conditions. Aggregates of clay were defined by geometrical analysis of face-to-face proximity of platelets with size distribution described by a log-normal function. The isotropy of the microstructure was assessed by computing a scalar order parameter. The number of platelets per aggregate and anisotropy of the microstructure both increases with platelet plan area. The system becomes more ordered and aggregate size increases with increasing pressure until maximum ordered state at confining pressure of 50 atm. Further increase of pressure slides platelets relative to each other leading to smaller aggregate size. The results show aggregate size of (3–8) platelets for sodium-smectite in agreement with experiments (3–10). The geometrical arrangement of aggregates affects mechanical properties of the system. The elastic properties of the meso-scale aggregate assembly are reported and compared with nanoindentation experiments. It is found that the elastic properties at this scale are close to the cubic systems. The elastic stiffness and anisotropy of the assembly increases with the size of the platelets and the level of external pressure.

Structure-Based Multiscale Approach for Identification of Interaction Partners of PDZ Domains

PDZ domains are peptide recognition modules which mediate specific protein-protein interactions and are known to have a complex specificity landscape. We have developed a novel structure-based multiscale approach which identifies crucial specificity determining residues (SDRs) of PDZ domains from explicit solvent molecular dynamics (MD) simulations on PDZ-peptide complexes and uses these SDRs in combination with knowledge-based scoring functions for proteomewide identification of their interaction partners. Multiple explicit solvent simulations ranging from 5 to 50 ns duration have been carried out on 28 PDZ-peptide complexes with known binding affinities. MM/PBSA binding energy values calculated from these simulations show a correlation coefficient of 0.755 with the experimental binding affinities. On the basis of the SDRs of PDZ domains identified by MD simulations, we have developed a simple scoring scheme for evaluating binding energies for PDZ-peptide complexes using residue based statistical pair potentials. This multiscale approach has been benchmarked on a mouse PDZ proteome array data set by calculating the binding energies for 217 different substrate peptides in binding pockets of 64 different mouse PDZ domains. Receiver operating characteristic (ROC) curve analysis indicates that, the area under curve (AUC) values for binder vs nonbinder classification by our structure based method is 0.780. Our structure based method does not require experimental PDZ-peptide binding data for training.

Specific Non-Local Interactions Are Not Necessary for Recovering Native Protein Dynamics

The elastic network model (ENM) is a widely used method to study native protein dynamics by normal mode analysis (NMA). In ENM we need information about all pairwise distances, and the distance between contacting atoms is restrained to the native value. Therefore ENM requires O(N2) information to realize its dynamics for a protein consisting of N amino acid residues. To see if (or to what extent) such a large amount of specific structural information is required to realize native protein dynamics, here we introduce a novel model based on only O(N) restraints. This model, named the ‘contact number diffusion’ model (CND), includes specific distance restraints for only local (along the amino acid sequence) atom pairs, and semi-specific non-local restraints imposed on each atom, rather than atom pairs. The semi-specific non-local restraints are defined in terms of the non-local contact numbers of atoms. The CND model exhibits the dynamic characteristics comparable to ENM and more correlated with the explicit-solvent molecular dynamics simulation than ENM. Moreover, unrealistic surface fluctuations often observed in ENM were suppressed in CND. On the other hand, in some ligand-bound structures CND showed larger fluctuations of buried protein atoms interacting with the ligand compared to ENM. In addition, fluctuations from CND and ENM show comparable correlations with the experimental B-factor. Although there are some indications of the importance of some specific non-local interactions, the semi-specific non-local interactions are mostly sufficient for reproducing the native protein dynamics.

Transition Paths of Met-Enkephalin from Markov State Modeling of a Molecular Dynamics Trajectory

Conformational states and their interconversion pathways of the zwitterionic form of the pentapeptide Met-enkephalin (MetEnk) are identified. An explicit solvent molecular dynamics (MD) trajectory is used to construct a Markov state model (MSM) based on dihedral space clustering of the trajectory, and transition path theory (TPT) is applied to identify pathways between open and closed conformers. In the MD trajectory, only four of the eight backbone dihedrals exhibit bistable behavior. Defining a conformer as the string XXXX with X = "+" or "-" denoting, respectively, positive or negative values of a given dihedral angle and obtaining the populations of these conformers shows that only four conformers are highly populated, implying a strong correlation among these dihedrals. Clustering in dihedral space to construct the MSM finds the same four bistable dihedral angles. These state populations are very similar to those found directly from the MD trajectory. TPT is used to obtain pathways, parametrized by committor values, in dihedral state space that are followed in transitioning from closed to open states. Pathway costs are estimated by introducing a kinetics-based procedure that orders pathways from least (shortest) to greater cost paths. The least costly pathways in dihedral space are found to only involve the same XXXX set of dihedral angles, and the conformers accessed in the closed to open transition pathways are identified. For these major pathways, a correlation between reaction path progress (committors) and the end-to-end distance is identified. A dihedral space principal component analysis of the MD trajectory shows that the first three modes capture most of the overall fluctuation, and pick out the same four dihedrals having essentially all the weight in those modes. A MSM based on root-mean-square backbone clustering was also carried out, with good agreement found with dihedral clustering for the static information, but with results that differ significantly for the pathway analysis.

Structural and Dynamic Properties of the Human Prion Protein

Prion diseases involve the conformational conversion of the cellular prion protein (PrPC) to its misfolded pathogenic form (PrPSc). To better understand the structural mechanism of this conversion, we performed extensive all-atom, explicit-solvent molecular-dynamics simulations for three structures of the wild-type human PrP (huPrP) at different pH values and temperatures. Residue 129 is polymorphic, being either Met or Val. Two of the three structures have Met in position 129 and the other has Val. Lowering the pH or raising the temperature induced large conformational changes of the C-terminal globular domain and increased exposure of its hydrophobic core. In some simulations, HA and its preceding S1-HA loop underwent large displacements. The C-terminus of HB was unstable and sometimes partially unfolded. Two hydrophobic residues, Phe-198 and Met-134, frequently became exposed to solvent. These conformational changes became more dramatic at lower pH or higher temperature. Furthermore, Tyr-169 and the S2-HB loop, or the X-loop, were different in the starting structures but converged to common conformations in the simulations for the Met-129, but not the Val-129, protein. α-Strands and β-strands formed in the initially unstructured N-terminus. α-Strand propensity in the N-terminus was different between the Met-129 and Val129 proteins, but β-strand propensity was similar. This study reveals detailed structural and dynamic properties of huPrP, providing insight into the mechanism of the conversion of PrPC to PrPSc.

Atomistic mechanisms of huntingtin N‐terminal fragment insertion on a phospholipid bilayer revealed by molecular dynamics simulations

The huntingtin protein is characterized by a segment of consecutive glutamines (Q(N)) that is responsible for its fibrillation. As with other amyloid proteins, misfolding of huntingtin is related to Huntington's disease through pathways that can involve interactions with phospholipid membranes. Experimental results suggest that the N-terminal 17-amino-acid sequence (htt(NT)) positioned just before the Q(N) region is important for the binding of huntingtin to membranes. Through all-atom explicit solvent molecular dynamics simulations, we unveil the structure and dynamics of the htt(NT)Q(N) fragment on a phospholipid membrane at the atomic level. We observe that the insertion dynamics of this peptide can be described by four main steps-approach, reorganization, anchoring, and insertion-that are very diverse at the atomic level. On the membrane, the htt(NT) peptide forms a stable α-helix essentially parallel to the membrane with its nonpolar side-chains-mainly Leu-4, Leu-7, Phe-11 and Leu-14-positioned in the hydrophobic core of the membrane. Salt-bridges involving Glu-5, Glu-12, Lys-6, and Lys-15, as well as hydrogen bonds involving Thr-3 and Ser-13 with the phospholipids also stabilize the structure and orientation of the htt(NT) peptide. These observations do not significantly change upon adding the Q(N) region whose role is rather to provide, through its hydrogen bonds with the phospholipids' head group, a stable scaffold facilitating the partitioning of the htt(NT) region in the membrane. Moreover, by staying accessible to the solvent, the amyloidogenic Q(N) region could also play a key role for the oligomerization of htt(NT)Q(N) on phospholipid membranes.