Learning probabilistic models of hydrogen bond stability from molecular dynamics simulation trajectories

Background

The estimation of individual ancestry from genetic data has become essential to applied population genetics and genetic epidemiology. Software programs for calculating ancestry estimates have become essential tools in the geneticist's analytic arsenal.

Results

Here we describe four enhancements to ADMIXTURE, a high-performance tool for estimating individual ancestries and population allele frequencies from SNP (single nucleotide polymorphism) data. First, ADMIXTURE can be used to estimate the number of underlying populations through cross-validation. Second, individuals of known ancestry can be exploited in supervised learning to yield more precise ancestry estimates. Third, by penalizing small admixture coefficients for each individual, one can encourage model parsimony, often yielding more interpretable results for small datasets or datasets with large numbers of ancestral populations. Finally, by exploiting multiple processors, large datasets can be analyzed even more rapidly.

Conclusions

The enhancements we have described make ADMIXTURE a more accurate, efficient, and versatile tool for ancestry estimation.     

Learning generative models of molecular dynamics

We introduce three algorithms for learning generative models of molecular structures from molecular dynamics simulations. The first algorithm learns a Bayesian-optimal undirected probabilistic model over user-specified covariates (e.g., fluctuations, distances, angles, etc). L1 regularization is used to ensure sparse models and thus reduce the risk of over-fitting the data. The topology of the resulting model reveals important couplings between different parts of the protein, thus aiding in the analysis of molecular motions. The generative nature of the model makes it well-suited to making predictions about the global effects of local structural changes (e.g., the binding of an allosteric regulator). Additionally, the model can be used to sample new conformations. The second algorithm learns a time-varying graphical model where the topology and parameters change smoothly along the trajectory, revealing the conformational sub-states. The last algorithm learns a Markov Chain over undirected graphical models which can be used to study and simulate kinetics. We demonstrate our algorithms on multiple molecular dynamics trajectories.     

Computer graphics presentations and analysis of hydrogen bonds from molecular dynamics simulation

To obtain better insights into the dynamic nature of hydrogen bonding, computer graphics representations were introduced as an aid for the analysis of molecular dynamics trajectories. A schematic representation of hydrogen bonding patterns is generated to reflect the frequency and the type of hydrogen bonding occurring during the simulation period. Various trajectory plots for monitoring geometrical parameters and for analyzing three-center hydrogen bonding were also generated. The methods proposed are applicable to a variety of biopolymers. In this study, hydrogen bonding in the d(G) · d(C)₆ system was examined. For the nucleic acid fragments examined, three-center hydrogen bonds can be classified as in-plane and major or minor groove types. The in-plane three-center hydrogen bond represents a stable state in which both bonds simultaneously satisfy the relaxed hydrogen bonding criteria for a measurable period. On the other hand, groove three-center hydrogen bonds behave as a transient intermediate state in a flip-flop hydrogen bonding system.     

Fuzzy cluster analysis of molecular dynamics trajectories.

We propose fuzzy clustering as a method to analyze molecular dynamics (MD) trajectories, especially of proteins and polypeptides. A fuzzy cluster analysis locates classes of similar three-dimensional conformations explored during a molecular dynamics simulation. The method can be readily applied to results from both equilibrium and nonequilibrium simulations, with clustering on either global or local structural parameters. The potential of this technique is illustrated by results from fuzzy cluster analyses of trajectories from MD simulations of various fragments of human parathyroid hormone (PTH). For large molecules, it is more efficient to analyze the clustering of root-mean-square distances between conformations comprising the trajectory. We found that the results of the clustering analysis were unambiguous, in terms of the optimal number of clusters of conformations, for the majority of the trajectories examined. The conformation closest to the cluster center can be chosen as being representative of the class of structures making up the cluster, and can be further analyzed, for example, in terms of its secondary structure. The CPU time used by the cluster analysis was negligible compared to the MD simulation time.     

Real-time interactive frequency filtering of molecular dynamics trajectories

Molecular dynamics simulations of atomic motion in protein and nucleic acid molecules must be done on a femtosecond time-scale. Much of this rapid motion is unimportant for the slower changes that are most relevant to biological function (conformational changes, substrate binding, protein folding). The high-frequency motion makes simulations computationally expensive. More importantly, the high frequencies obscure visualization of the relevant dynamics processes. Sessions, Dauber-Osguthorpe and Osguthorpe presented a method for removing high-frequency motions from atomic co-ordinates of trajectories generated by simulation. While that study used fast Fourier methods and emphasized the use of filtering for analysis of trajectories, this communication describes a new method that makes it much easier to use frequency filtering in programs that display trajectories as a sequence of moving images. Tests of the method on systems extending from pure water to proteins and nucleic acid molecules in vacuo and in solution have demonstrated its general utility. Impressed with the power and simplicity of the new method, we wish to present it in sufficient detail to allow others to implement it themselves.     

Implicit solvent models for flexible protein-protein docking by molecular dynamics simulation

The suitability of three implicit solvent models for flexible protein-protein docking by procedures using molecular dynamics simulation is investigated. The three models are (i) the generalized Born (GB) model implemented in the program AMBER6.0; (ii) a distance-dependent dielectric (DDD) model; and (iii) a surface area-dependent model that we have parameterized and call the NPSA model. This is a distance-dependent dielectric model modified by neutralizing the ionizable side-chains and adding a surface area-dependent solvation term. These solvent models were first tested in molecular dynamics simulations at 300 K of the native structures of barnase, barstar, segment B1 of protein G, and three WW domains. These protein structures display a range of secondary structure contents and stabilities. Then, to investigate the performance of the implicit solvent models in protein docking, molecular dynamics simulations of barnase/barstar complexation, as well as PIN1 WW domain/peptide complexation, were conducted, starting from separated unbound structures. The simulations show that the NPSA model has significant advantages over the DDD and GB models in maintaining the native structures of the proteins and providing more accurate docked complexes.     

Structure and stability of the insulin dimer investigated by molecular dynamics simulation

Molecular dynamics simulation indicates that the dynamical behaviour of the insulin dimer is asymmetric. Atomic level knowledge of the interaction modes and protein conformation in the solvation state identifies dynamical structures, held by hydrogen bonds that stabilize, mainly in one monomer, the interaction between the chains. Dynamic cross-correlation analysis shows that the two insulin monomers behave asymmetrically and are almost independent. Solvation energy, calculated to evaluate the contribute of each interface residue to the dimer association pattern, well compares with the experimental association state found in protein mutants indicating that this parameter is an important factor to explain the association properties of mutated insulin dimers.     

Effect of hydrogen coverage on the buckling of penta-graphene by molecular dynamics simulation

By employing a series of MD simulations, buckling behaviour of penta-graphene and functionalised penta-graphene having different hydrogen (H) coverage is presented in this study. To this end, the buckling onset strain is determined for different systems. The results reveal that the new allotrope is slightly stiffer than graphene. Moreover, the effect of H adatoms in the range 0–56% on buckling behaviour is investigated. Finally it is shown that the H coverage has an influence on stability, and ripple-type destortion of the penta-graphene under compression.     

Data mining of molecular dynamics trajectories of nucleic acids

Analysis, storage, and transfer of molecular dynamic trajectories are becoming the bottleneck of computer simulations. In this paper we discuss different approaches for data mining and data processing of huge trajectory files generated from molecular dynamic simulations of nucleic acids.     

Sequence specific hydrogen bond of DNA with denaturants affects its stability: Spectroscopic and simulation studies

BackgroundSeveral different small molecules have been used to target the DNA helix in order to treat the diseases caused by its mutation. Guanidinium(Gdm⁺) and urea based drugs have been used for the diseases related to central nervous system, also as the anti-inflammatory and chemotherapeutic agent. However, the role of Gdm⁺ and urea in the stabilization/destabilization of DNA is not well understood.MethodsSpectroscopic techniques along with molecular dynamics (MD) simulation have been performed on different sequences of DNA in the presence of guanidinium chloride (GdmCl) and urea to decode the binding of denaturants with DNA and the role of hydrogen bond with the different regions of DNA in its stability/destability.Results and conclusionOur study reveals that, Gdm⁺ of GdmCl and urea both intrudes into the groove region of DNA along with the interaction with its phosphate backbone. However, interaction of Gdm⁺ and urea with the nucleobases in the groove region is different. Gdm⁺ forms the intra-strand hydrogen bond with the central region of the both sequences of DNA whereas inter-strand hydrogen bond along with water assisted hydrogen bond takes place in the case of urea. The intra-strand hydrogen bond formation capability of Gdm⁺ with the nucleobases in the minor groove of DNA decreases its groove width which probably causes the stabilization of B-DNA in GdmCl. In contrast, the propensity of the formation of inter-strand hydrogen bond of urea with the nucleobases in the groove region of DNA without affecting the groove width destabilizes B-DNA as compared to GdmCl. This study depicts that the opposite effect of GdmCl and urea on the stability is a general property of B-DNA. However, the extent of stabilization/destabilization of DNA in Gdm⁺ and urea depend on its sequence probably due to the difference in the intra/inter-strand hydrogen bonding with different bases present in both the sequences of DNA.General significanceThe information obtained from this study will be useful for the designing of Gdm⁺ based drug molecule which can target the DNA more specifically and selectively.     

The asparagine-stabilized beta-turn of apamin: contribution to structural stability from dynamics simulation and amide hydrogen exchange analysis

Molecular dynamics simulations of bee venom apamin, and an analogue having an Asn to Ala substitution at residue 2 (apamin-N2A), were analyzed to explore the contribution of hydrogen bonds involving Asn2 to local (beta-turn residues N2, C3, K4, A5) and global stability. The wild-type peptide retained a stable conformation during 2.4 ns of simulation at 67 degrees C, with high beta-turn stability characterized by backbone-side chain hydrogen bonds involving beta-turn residues K4 and A5, with the N2 side chain amide carbonyl. The loss of stabilizing interactions involving the N2 side chain resulted in the loss of the beta-turn conformation in the apamin N2A simulations (27 or 67 degrees C). This loss of beta-turn stability propagates throughout the peptide structure, with destabilization of the C-terminal helix connected to the N-terminal region by two disulfide bonds. Backbone stability in a synthetic peptide analogue (apamin-N2A) was characterized by NMR and amide hydrogen exchange measurements. Consistent with the simulations, loss of hydrogen bonds involving the N2 side chain resulted in destabilization of both the N-terminal beta-turn and the C-terminal helix. Amide exchange protection factors in the C-terminal helix were reduced by 9-11-fold in apamin N2A as compared with apamin, corresponding to free energy (deltaDeltaG(uf)) of around 1.5 kcal M(-1) at 20 degrees C. This is equivalent to the contribution of hydrogen bond interactions involving the N2 side chain to the stability of the beta-turn. Together with additional measures of exchange protection factors, the three main contributions to backbone stability in apamin that account for virtually the full thermodynamic stability of the peptide have been quantitated.     

Importance of explicit salt ions for protein stability in molecular dynamics simulation.

The accurate and efficient treatment of electrostatic interactions is one of the challenging problems of molecular dynamics simulation. Truncation procedures such as switching or shifting energies or forces lead to artifacts and significantly reduced accuracy. The particle mesh Ewald (PME) method is one approach to overcome these problems by providing a computationally efficient means of calculating all long-range electrostatic interactions in a periodic simulation box by use of fast Fourier transformation techniques. For the application of the PME method to the simulation of a protein with a net charge in aqueous solution, counterions are added to neutralize the system. The usual procedure is to add charge-balancing counterions close to charged residues to neutralize the protein surface. In the present article, we show that for MD simulation of a small protein of marginal stability, the YAP-WW domain, explicit modeling of 0.2 M ionic strength (in addition to the charge-balancing counterions) is necessary to maintain a stable protein structure. Without explicit ions throughout the periodic simulation box, the charge-balancing counterions on the protein surface diffuse away from the protein, resulting in destruction of the beta-sheet secondary structure of the WW domain.     

Structural stability of myoglobin and glycomyoglobin: a comparative molecular dynamics simulation study

Glycoproteins are formed as the result of enzymatic glycosylation or chemical glycation in the body, and produced in vitro in industrial processes. The covalently attached carbohydrate molecule(s) confer new properties to the protein, including modified stability. In the present study, the structural stability of a glycoprotein form of myoglobin, bearing a glucose unit in the N-terminus, has been compared with its native form by the use of molecular dynamics simulation. Both structures were subjected to temperatures of 300 and 500 K in an aqueous environment for 10 ns. Changes in secondary structures and RMSD were then assessed. An overall higher stability was detected for glycomyoglobin, for which the most stable segments/residues were highlighted and compared with the native form. The simple addition of a covalently bound glucose is suggested to exert its stabilizing effect via increased contacts with surrounding water molecules, as well as a different pattern of interactions with neighbor residues.

Electronic supplementary material

The online version of this article (doi:10.1007/s10867-015-9383-2) contains supplementary material, which is available to authorized users.     

Conformational stability of PCID2 upon DSS1 binding with molecular dynamics simulation

DSS1 is a small acidic intrinsically disordered protein (IDP) that can fold upon binding with PCID2 TREX-2. The resulting complex plays a key role in mRNA export. However, the binding mechanism between DSS1 and PCID2 is unsolved. Here, three independent 500-ns molecular dynamics (MD) simulations were performed to study the DSS1–PCID2 binding mechanism by comparing apo-PCID2 and bound PCID2. The results show that the conformational variation of bound PCID2 is smaller than that of apo-PCID2, especially in the binding domain of two helices (helix IV and VIII). The probability of coil formation between helix III and helix IV of bound PCID2 increases, and a short anti-parallel β-sheet forms upon DSS1 binding. The decomposition of binding free energy into protein and residue pairs suggests that electrostatic and hydrophobic interactions play key roles in the recognition between DSS1 and PCID2. There is a hydrophobic core of seven residues in DSS1 favorable to the binding of PCID2. These analytical methods can be used to reveal the recognition mechanisms of other IDPs and their partners.

Open image in new window

Figure Differences of second structure of PCID2 in bound and unbound states. The interaction surface between the helix VIII of PCID2 and helix of DSS1

     

Congruent qualitative behavior of complete and reconstructed phase space trajectories from biomolecular dynamics simulation

Central to the study of a complex dynamical system is knowledge of its phase space behavior. Experimentally, it is rarely possible to record a system's (multidimensional) phase space variables. Rather, the system is observed via one (or few) scalar-valued signal(s) of emission or response. In dynamical systems analysis, the multidimensional phase space of a system can be reconstructed by manipulation of a one-dimensional signal. The trick is in the construction of a (higher-dimensional) space through the use of a time lag (or delay) on the signal time series. The trajectory in this embedding space can then be examined using phase portraits generated in selected subspaces. By contrast, in computer simulation, one has an embarrassment of riches: direct access to the complete multidimensional phase space variables, at arbitrary time resolution and precision. Here, the problem is one of reducing the dimensionality to make analysis tractable. This can be achieved through linear or nonlinear projection of the trajectory into subspaces containing high information content. This study considers trajectories of the small protein crambin from molecular dynamics simulations. The phase space behavior is examined using principal component analysis on the Cartesian coordinate covariance matrix of 138 dimensions. In addition, the phase space is reconstructed from a one dimensional signal, representing the radius of gyration of the structure along the trajectory. Comparison of low-dimensional phase portraits obtained from the two methods shows that the complete phase space distribution is well represented by the reconstruction. The study suggests that it may be possible to develop a deeper connection between the experimental and simulated dynamics of biomolecules via phase space reconstruction using data emerging from recent advances in single-molecule time-resolved biophysical techniques.     

Insights into structure, stability, and toxicity of monomeric and aggregated polyglutamine models from molecular dynamics simulations

Nine genetically inherited neurodegenerative diseases are linked to abnormal expansions of a polyglutamine (polyQ) encoding region. Over the years, several structural models for polyQ regions have been proposed and confuted. The cross-β-spine steric zipper motif, identified recently for the GNNQQNY peptide, represents an attractive model for amyloid fibers formed by polyQ fragments. Here we report a detailed molecular dynamics investigation of polyQ models assembled by cross-β-spine steric zipper motifs. Our simulations indicate clearly that these assemblies are very stable. Glutamine side chains contribute strongly to the overall stability of the models by fitting perfectly within the zipper. In contrast to GNNQQNY zipper motifs, hydrogen bonding interactions provide a significant contribution to the overall stability of polyQ models. Molecular dynamics simulations carried out on monomeric polyQ forms (composed by 40-60 residues) show clearly that they can also assume structures stabilized by steric zipper motifs. Based on these findings, we build monomeric polyQ models that can explain recent data on the toxicity exerted by these species. In a more general context, our data suggests that polyQ models with interdigitated side chains can provide a structural rationale to several literature experiments on polyQ formation, stability, and toxicity.     

Simulation of nitroxide electron paramagnetic resonance spectra from brownian trajectories and molecular dynamics simulations

A simulated continuous wave electron paramagnetic resonance spectrum of a nitroxide spin label can be obtained from the Fourier transform of a free induction decay. It has been previously shown that the free induction decay can be calculated by solving the time-dependent stochastic Liouville equation for a set of Brownian trajectories defining the rotational dynamics of the label. In this work, a quaternion-based Monte Carlo algorithm has been developed to generate Brownian trajectories describing the global rotational diffusion of a spin-labeled protein. Also, molecular dynamics simulations of two spin-labeled mutants of T4 lysozyme, T4L F153R1, and T4L K65R1 have been used to generate trajectories describing the internal dynamics of the protein and the local dynamics of the spin-label side chain. Trajectories from the molecular dynamics simulations combined with trajectories describing the global rotational diffusion of the protein are used to account for all of the dynamics of a spin-labeled protein. Spectra calculated from these combined trajectories correspond well to the experimental spectra for the buried site T4L F153R1 and the helix surface site T4L K65R1. This work provides a framework to further explore the modeling of the dynamics of the spin-label side chain in the wide variety of labeling environments encountered in site-directed spin labeling studies.     

The self-diffusion and hydrogen bond interaction in neat liquid alkanols: a molecular dynamic simulation study

Self-diffusion of methanol, ethanol, 1-propanol and 2-propanol has been studied by molecular dynamics simulation in the temperature range between the melting pressure curve and 478 K at pressures up to 300 MPa. The simulation results on self-diffusion of methanol, ethanol and 2-propanol (for 2-propanol, at high temperatures) agree well with experiment, which suggests that the simulation method is a powerful tool to obtain self-diffusion coefficients over wide range of temperature and pressure, under which it is rather difficult for experiments. The local structures of methanol, ethanol and 2-propanol are investigated by calculating the radial distribution functions, H-bond numbers, coordination numbers and the ratios of H-bond number divided by coordination number. The correlation between self-diffusion and structural properties, and the influence of temperature and pressure on them are discussed. The degree of forming H-bond space network in methanol, ethanol and water is higher than that in 2-propanol, and they are all higher than those in ammonia and methylamine. The simulation results demonstrate that the effect of hydrogen bonding on the translational dynamics in methanol and ethanol is more pronounced than that in 2-propanol.