Search

The prokaryotic monotopic membrane protein squalene-hopene cyclase (SHC) is homologous to a human enzyme responsible for cholesterol formation. Using molecular dynamics in explicit water, a single monomer of SHC was simulated using the GROMOS 45A3 force field, once in complex with an inhibitor and once in an uncomplexed form. The protein exhibits significant stability on the level of secondary and tertiary structure even outside of its native membrane environment. Analysis of the fluctuations of the complexed and the uncomplexed SHC confirms the previously made suggestions for the ligand entrance channel and reveals some of its novel dynamical features that might be of functional importance. To examine the potential of computationally designing SHC ligands and study their thermodynamics of binding, the relative free energies of binding of a series of structurally similar anticholesteremic inhibitors of SHC were calculated using single-step perturbation (SSP) and thermodynamic integration (TI) techniques. While neither technique succeeds in quantitatively matching the relatively small experimental values, TI qualitatively reproduces the relative order of the experimental affinities, but SSP does not. Detailed comparisons and potential reasons for this are given.

The common characteristics of amyloid and amyloid-like fibrils from disease- and non-disease-associated proteins offer the prospect that well-defined model systems can be used to systematically dissect the driving forces of amyloid formation. We recently reported the de novo designed cc beta peptide model system that forms a native-like coiled-coil structure at low temperatures and which can be switched to amyloid-like fibrils by increasing the temperature. Here, we report a detailed molecular description of the system in its fibrillar state by characterizing the cc beta-Met variant using several microscopic techniques, circular dichroism spectroscopy, X-ray fiber diffraction, solid-state nuclear magnetic resonance, and molecular dynamics calculations. We show that cc beta-Met forms amyloid-like fibrils of different morphologies on both the macroscopic and atomic levels, which can be controlled by variations of assembly conditions. Interestingly, heterogeneity is also observed along single fibrils. We propose atomic models of the cc beta-Met amyloid-like fibril, which are in good agreement with all experimental data. The models provide a rational explanation why oxidation of methionine residues completely abolishes cc beta-Met amyloid fibril formation, indicating that a small number of site-specific hydrophobic interactions can play a major role in the packing of polypeptide-chain segments within amyloid fibrils. The detailed structural information available for the cc beta model system provides a strong molecular basis for understanding the influence and relative contribution of hydrophobic interactions on native-state stability, kinetics of fibril formation, fibril packing, and polymorphism.

Over the last 30 years, computation based on molecular models is playing an increasingly important role in biology, biological chemistry and biophysics. Since only a very limited number of properties of biomolecular systems are actually accessible to measurement by experimental means, computer simulation complements experiments by providing not only averages, but also distributions and time series of any definable, observable or non-observable, quantity. Biomolecular simulation may be used (i) to interpret experimental data, (ii) to provoke new experiments, (iii) to replace experiments and (iv) to protect intellectual property. Progress over the last 30 years is sketched and perspectives are outlined for the future.

Methods to search for low-energy conformations, to generate a Boltzmann-weighted ensemble of configurations, or to generate classical-dynamical trajectories for molecular systems in the condensed liquid phase are briefly reviewed with an eye to application to biomolecular systems. After having chosen the degrees of freedom and method to generate molecular configurations, the efficiency of the search or sampling can be enhanced in various ways: (i) efficient calculation of the energy function and forces, (ii) application of a plethora of search enhancement techniques, (iii) use of a biasing potential energy term, and (iv) guiding the sampling using a reaction or transition pathway. The overview of the available methods should help the reader to choose the combination that is most suitable for the biomolecular system, degrees of freedom, interaction function, and molecular or thermodynamic properties of interest., ((c) 2007 Wiley Periodicals, Inc. J Comput Chem, 2008.)

Well-established semi-empirical quantum-mechanical methods such as AM1, PM3 and MNDO poorly describe hydrogen-bonding interactions. The recently developed OM1, OM2 and OM3 methods, which include explicit corrections in the one-electron terms of the Hamiltonian to counteract for the neglect of orthogonalization effects, show an improved description of hydrogen bonding. In the current work we compare the performance of the OMx models on the one hand and the AM1, PM3 and MNDO methods (as well as the SCC-DFTB method) on the other when used in combination with a classical force field in QM/MM calculations on a water dimer in the gas phase, and in QM/MM simulations of a QM water 'solute' in liquid MM water. It is found that when using the OMx methods to describe the QM water molecule, values for QM-MM interaction energies and the molecular dipole moment of the QM solute are significantly closer to experiment than when using the other semi-empirical methods. In addition, the compatibility of OM3 with the classical Hamiltonian improves upon explicitly including electronic polarization in the MM subsystem.

The temperature dependence of thermodynamic quantities, such as heat capacity, entropy and free enthalpy, may be obtained by using well-known equations that relate these quantities to the enthalpy of the molecular system of interest at a range of temperatures. In turn, the enthalpy of a molecular system can be estimated from molecular dynamics simulations of an appropriate model. To demonstrate this, we have investigated the temperature dependence of the enthalpy, heat capacity, entropy and free enthalpy of a system that consists of a beta-heptapeptide in methanol and have used the statistical mechanics relationships to describe the thermodynamics of the folding/unfolding equilibrium of the peptide. The results illustrate the power of current molecular simulation force fields and techniques in establishing the link between thermodynamic quantities and conformational distributions.

Introducing experimental values as restraints into molecular dynamics (MD) simulation to bias the values of particular molecular properties, such as nuclear Overhauser effect intensities or distances, dipolar couplings, 3 J-coupling constants, chemical shifts or crystallographic structure factors, towards experimental values is a widely used structure refinement method. Because multiple torsion angle values varphi correspond to the same 3 J-coupling constant and high-energy barriers are separating those, restraining 3 J-coupling constants remains difficult. A method to adaptively enforce restraints using a local elevation (LE) potential energy function is presented and applied to 3 J-coupling constant restraining in an MD simulation of hen egg-white lysozyme (HEWL). The method succesfully enhances sampling of the restrained torsion angles until the 37 experimental 3 J-coupling constant values are reached, thereby also improving the agreement with the 1,630 experimental NOE atom-atom distance upper bounds. Afterwards the torsional angles varphi are kept restrained by the built-up local-elevation potential energies.

Widely used programs for molecular dynamics simulation of (bio)molecular systems are the Verlet and leapfrog algorithms. In these algorithms, the particle velocities are less accurately propagated than the positions. Important quantities for the simulation such as the temperature and the pressure involve the squared velocities at full time steps. Here, we derive an expression for the squared particle velocity at full time step in the leapfrog scheme, which is more accurate than the standardly used one. In particular, this allows us to show that the full time step kinetic energy of a particle is more accurately computed as the average of the kinetic energies at previous and following half steps than as the square of the average velocity as implemented in various molecular dynamics codes. Use of the square of the average velocity introduces a systematic bias in the calculation of the instantaneous temperature and pressure of a molecular dynamics system. We show the consequences when the system is coupled to a thermostat and a barostat.

An expression for the atomic forces in simulations using the charge-on-spring (COS) polarizable model is rederived. In previous implementations of COS-based models, contributions arising from the dependence of the induced dipoles (i.e., the positions of the charges-on-spring) on the coordinates of the other sites in the system were not taken into account. However, from calculations on gas-phase dimers we found a significant contribution of these terms. Errors in the forces when neglecting these contributions in condensed-phase calculations can be significantly reduced by choosing an appropriately large value for the size of the charge-on-spring.

A comparison of different treatments of bond-stretching interactions in molecular dynamics simulation is presented. Relative free energies from simulations using rigid bonds maintained with the SHAKE algorithm, using partially rigid bonds maintained with a recently introduced flexible constraints algorithm, and using fully flexible bonds are compared in a multi-configurational thermodynamic integration calculation of changing liquid water into liquid methanol. The formula for the free energy change due to a changing flexible constraint in a flexible constraint simulation is derived. To allow for a more direct comparison between these three methods, three different pairs of models for water and methanol were used: a flexible model (simulated without constraints and with flexible constraints), a rigid model (simulated with standard hard constraints), and an alternative flexible model (simulated with flexible constraints and standard hard constraints) in which the ideal or constrained bond lengths correspond to the average bond lengths obtained from a short simulation of the unconstrained flexible model. The particular treatment of the bonds induces differences of up to 2 % in the liquid densities, whereas (excess) free energy differences of up to 5.7 (4.3) kJ mol(-1) are observed. These values are smaller than the differences observed between the three different pairs of methanol/water models: up to 5 % in density and up to 8.5 kJ mol(-1) in (excess) free energy.

Schizophyllan is a beta(1-->3)-D-glucan polysaccharide with beta(1-->6)-branched lateral glucose residues that presents a very stiff triple-helical structure under most experimental conditions. Despite the remarkable stability of this structure (which persists up to 120 degrees C in aqueous solution), schizophyllan undergoes a major change of state around 7 degrees C in water that has been hypothesized to result from an order-disorder transition in the lateral residues. This hypothesis is only supported by indirect experimental evidence and detailed knowledge (at the atomic level) concerning hydrogen-bonding networks, interactions with the solvent molecules, orientational freedom of the lateral residues, and orientational correlations among them is still lacking. In this study explicit-solvent molecular dynamics simulations of a schizophyllan fragment (complemented by simulations of its tetrasaccharide monomer) are performed at three different temperatures (273 K, 350 K, and 450 K) and with two different types of boundary conditions (finite nonperiodic or infinite periodic fragment) as an attempt to provide detailed structural and dynamical information about the triple-helical conformation in solution and the mechanism of the low-temperature transition. These simulations suggest that three important driving forces for the high stability of the triple helix are i), the limited conformational work involved in its formation; ii), the formation of a dense hydrogen-bonding network at its center; and iii), the formation of interchain hydrogen bonds between main-chain and lateral glucose residues. However, these simulations evidence a moderate and continuous variation of the simulated observables upon increasing the temperature, rather than a sharp transition between the two lowest temperatures (that could be associated with the state transition). Although water-mediated hydrogen-bonded association of neighboring lateral residues is observed, this interaction is not strong enough to promote the formation of an ordered state (correlated motions of the lateral residues), even at the lowest temperature considered.

In a previous combined QM/MM molecular dynamics (MD) study from our laboratory on the identity SN2 reaction between a chloride anion and an amino chloride in liquid dimethyl ether (DME), an increase in the free energy activation barrier was observed in the condensed phase when compared to the gas-phase activation energy. Here we reproduce these findings, but when comparing the condensed-phase potential of mean force (PMF) with the free energy profile in the gas phase (obtained from Monte Carlo simulations), we observe a smaller solvent effect on the activation barrier of the reaction. In a next step, we introduce an explicit description of electronic polarization in the MM (solvent) part of the system. A polarizable force field for liquid DME was developed based on the charge-on-spring (COS) model, which was calibrated to reproduce thermodynamic properties of the nonpolarizable model in classical MD simulations. The COS model was implemented into the MNDO/GROMOS interface in a special version of the QM/MM software ChemShell, which was used to investigate the effect of solvent polarization on the free energy profile of the reaction under study. A higher activation barrier was obtained using the polarizable solvent model than with the nonpolarizable force field, due to a better solvation of and a stronger polarization of solvent molecules around the separate reactants. The obtained PMFs were subjected to an energy-entropy decomposition of the relative solvation free energies of the reactant complex along the reaction coordinate, to investigate in a quantitative manner whether the solvent (polarization) effects are mainly due to favorable QM-MM (energetic) interactions.

Detailed knowledge of how networks of surface salt bridges contribute to protein thermal stability is essential not only to understand protein structure and function but also to design thermostable proteins for industrial applications. Experimental studies investigating thermodynamic stability through measurements of free energy associated with mutational alterations in proteins provide only macroscopic evidence regarding the structure of salt-bridge networks and assessment of their contribution to protein stability. Using explicit-solvent molecular dynamics simulations to provide insight on the atomic scale, we investigate here the structural stability, defined in terms of root-mean-square fluctuations, of a short polypeptide designed to fold into a stable trimeric coiled coil with a well-packed hydrophobic core and an optimal number of intra- and interhelical surface salt bridges. We find that the increase of configurational entropy of the backbone and side-chain atoms and decreased pair correlations of these with increased temperature are consistent with nearly constant atom-positional root-mean-square fluctuations, increased salt-bridge occupancies, and stronger electrostatic interactions in the coiled coil. Thus, our study of the coiled coil suggests a mechanism in which well-designed salt-bridge networks could accommodate stochastically the disorder of increased thermal motion to produce thermostability.

The lack of an explicit description of electronic polarization in nonpolarizable force fields usually results in an incomplete transferability of force-field parameter sets when applied in simulations of the system of interest in either a polar or an apolar environment. For example, the use of nonpolarizable parameter sets optimized to reproduce experimental data on properties of pure liquids of polar compounds commonly yields too low solubilities in water for the corresponding compounds. The reason is that the fixed charge distributions calibrated for the pure liquid might correspond to too low molecular dipole moments in case of hydration. In the current study, we quantitatively show that explicit inclusion of electronic polarization can improve the transferability of biomolecular force-field parameter sets. With this aim, free energies of polarization, DeltaGpola, have been calculated, with DeltaGpola corresponding to the free energy difference between identical systems described by a polarizable and a nonpolarizable model. Using a nonpolarizable model and a polarizable one (based on the charge-on-spring approach) for dimethyl ether (DME), which were both parametrized to reproduce experimental values for pure liquid properties, small values were found for DeltaGpola for the pure liquid or when a DME solute was solvated in the apolar solvent cyclohexane. For the solute hydrated in water, however, DeltaGpola was found to be of the same order of magnitude as the discrepancy between the free energy of hydration from simulation using a nonpolarizable solute model and the experimental value. Thus, introducing polarizabilities clearly improves the transferability of the parameter set. Additionally, in calculations of an anion solvated in DME, DeltaGpola for the solvent adopted relatively large values. From an estimation of the errors in the calculated free energy differences, it was furthermore shown that the calculation of DeltaGpola offers an effective and accurate method to obtain differences in solvation (or excess) free energies between systems described by polarizable and nonpolarizable models when compared to a direct calculation of solvation (or excess) free energies.

The polypyrimidine tract binding protein (PTB) is involved in many physiological processes, including alternative splicing, internal ribosomal entry side (IRES)-mediated initiation of translation, and polyadenylation, as well as in ensuring mRNA stability. However, the role of PTB in these processes is not fully understood, and this has motivated us to undertake a computational study of the protein. PTB RNA binding domains (RBDs) 3 and 4 and their complexes with oligopyrimidine RNAs were simulated using the GROMOS simulation software using the GROMOS 45A4 force field. First, the stability and fluctuations of the tertiary fold and of the secondary structural elements in individual domains, the combined RBD34 domain, and their complexes with RNA were studied. Second, the simulation results were validated against the experimental NMR NOE data. The analysis of hydrogen bonding patterns, salt bridge networks, and stacking interactions of the RNA to the binding pockets of the protein domains showed that binding is not sequence-specific and that many RNA fragments can bind to them successfully. Further calculations of the relative free energy of binding for different polypyrimidine sequences were carried out using the thermodynamic integration (TI) and single-step perturbation (SSP) methods. It is was not possible to calculate the relative free energies with high accuracy, but the obtained results do give qualitative insights into PTB's affinity for different RNA sequences. Furthermore, the low-energy conformations of the complexes that were found provided additional information about the mechanism of binding.

A series of explicit-solvent molecular dynamics simulations of the protein ubiquitin are reported, which investigate the effect of environmental factors (presence of methanol cosolvent in the aqueous solution, neutral or low pH value, room or elevated temperature) on the structure, stability, and dynamics of the protein. The simulations are initiated either from the native structure of the protein or from a model of a partially folded state (A-state) that is known to exist at low pH in methanol-water mixtures. The main results of the simulations are: (1) The ubiquitin native structure is remarkably stable at neutral pH in water; (2) the addition of the methanol cosolvent enhances the stability of the secondary structure but weakens tertiary interactions within the protein; (3) this influence of methanol on the protein structure is enhanced at low pH, while the effect of lowering the pH in pure water is limited; and (4) the A-state of ubiquitin can be described as a set of relatively rigid secondary structure elements (a native-like beta-sheet and native-like alpha-helix plus two nonnative alpha-helices) connected by flexible linkers.

The authors present a method to calculate free energy differences between two states A and B "on the fly" from a single molecular dynamics simulation of a reference state R. No computer time has to be spent on the simulation of intermediate states. Only one state is sampled, i.e., the reference state R which is designed such that the subset of phase space important to it is the union of the parts of phase space important to A and B. Therefore, an accurate estimate of the relative free energy can be obtained by construction. The authors applied the method to four test systems (dipole inversion, van der Waals interaction perturbation, charge inversion, and water to methanol conversion) and compared the results to thermodynamic integration estimates. In two cases, the enveloping distribution sampling calculation was straightforward. However, in the charge inversion and the water to methanol conversion, Hamiltonian replica-exchange molecular dynamics of the reference state was necessary to observe transitions in the reference state simulation between the parts of phase space important to A and B, respectively. This can be explained by the total absence of phase space overlap of A and B in these two cases.

beta-Depsipeptides are beta-peptides in which one or more peptide linkages are replaced by ester linkages, resulting in a loss of a hydrogen-bond donor (N--H) and weakening of the corresponding carbonyl hydrogen-bond acceptor moiety. The effects of three of such peptide by ester substitutions in a hepta-beta-peptide upon its (un)folding equilibrium in methanol solution are investigated using molecular dynamics simulations and compared to experimental data from NMR spectroscopy. The simulated conformational ensembles largely reproduce the experimentally measured NOE and 3J-coupling constant data for the three different hepta-beta-peptides, and confirm the relative stabilities of the 3(14)-helical conformation, which is most weakened by substitution of the 4th peptide linkage and least by substitution of the 6th peptide linkage. The simulations are complementary to the experimental data by providing detailed insight into the conformational distributions that are compatible with the experimentally measured average values of observables., ((c) 2006 Wiley Periodicals, Inc.)

Thermodynamic data are often used to calibrate or test amomic-level (AL) force fields for molecular dynamics (MD) simulations. In contrast, the majority of coarse-grained (CG) force fields do not rely extensively on thermodynamic quantities. Recently, a CG force field for lipids, hydrocarbons, ions, and water, in which approximately four non-hydrogen atoms are mapped onto one interaction site, has been proposed and applied to study various aspects of lipid systems. To date, no extensive investigation of its capability to describe salvation thermodynamics has been undertaken. In the present study, a detailed picture of vaporization, solvation, and phase-partitioning thermodynamics for liquid hydrocarbons and water was obtained at CG and AL resolutions, in order to compare the two types or models and evaluate their ability to describe thermodynamic properties in the temperature range between 263 and 343 K. Both CG and AL models capture the experimental dependence of the thermodynamic properties on the temperature, albeit a systematically weaker dependence is found for the CG model. Moreover, deviations are found for solvation thermodynamics and for the corresponding enthalpy-entropy compensation for the CG model. Particularly water/oil repulsion seems to be overestimated. However, the results suggest that the thermodynamic properties considered should be reproducible by a CG model provided it is reparametrized on the basis of these liquid-phase properties.

Most processes occurring in a system are determined by the relative free energy between two or more states because the free energy is a measure of the probability of finding the system in a given state. When the two states of interest are connected by a pathway, usually called reaction coordinate, along which the free-energy profile is determined, this profile or potential of mean force (PMF) will also yield the relative free energy of the two states. Twelve different methods to compute a PMF are reviewed and compared, with regard to their precision, for a system consisting of a pair of methane molecules in aqueous solution. We analyze all combinations of the type of sampling (unbiased, umbrella-biased or constraint-biased), how to compute free energies (from density of states or force averaging) and the type of coordinate system (internal or Cartesian) used for the PMF degree of freedom. The method of choice is constraint-bias simulation combined with force averaging for either an internal or a Cartesian PMF degree of freedom.

The use of time-dependent restraints in molecular simulation in order to generate a conformational ensemble for molecules that is in accordance with measured ensemble averages for particular observable quantities is investigated. Using a model system consisting of liquid butane and the cyclic peptide antamanide the reproduction of particular average (3) J-coupling constant values in a molecular dynamics simulation is analysed. It is shown that the multiple-valuedness and the sizeable gradients of the Karplus curve relating (3) J-coupling constants measured in NMR experiments to the corresponding torsional-angle values cause severe problems when trying to restrain a (3) J-coupling constant to a value close to the extrema of the Karplus curve. The introduction of a factor oscillating with time into the restraining penalty function alleviates this problem and enhances the restrained conformational sampling.

Phosphodiesterases are a large class of enzymes mediating a number of physiological processes ranging from immune response to platelet aggregation to cardiac and smooth muscle relaxation. In particular, phosphodiesterase 5 (PDE5) plays an important role in mediating sexual arousal, and it is the central molecular target in treatments of erectile dysfunction. In this study, we look at the mechanism and thermodynamics of the binding of selective inhibitors sildenafil (Viagra) and vardenafil (Levitra) to PDE5 using molecular dynamics simulations. Our simulations of PDE5 with and without sildenafil suggest a binding mechanism in which two loops surrounding the binding pocket of the enzyme (H loop, residues 660-683, and M loop, 787-812) execute sizable conformational changes (∼1 nm), clamping the ligand in the pocket. Also, we note significant changes in the coordination pattern of the divalent ions in the active site of the enzyme, as well as marked changes in the collective motions of the enzyme when the ligand is bound. Using the thermodynamic integration approach we calculate the relative free energies of binding of sildenafil, vardenafil, and demethyl-vardenafil, providing a test of the quality of the force field and the ligand parametrization used. Finally, using the single-step perturbation (SSP) technique, we calculate the relative binding free energies of these three ligands as well. In particular, we focus on critical evaluation of the SSP technique and examine the effects of computational parallelization on the efficiency of the technique. As a technical improvement, we demonstrate that an ensemble of relatively short SSP trajectories (10 × 0.5 ns) markedly outperforms a single trajectory of the same total length (1 × 5 ns) when it comes to sampling efficiency, resulting in significant real-time savings.

Simulations of various beta-peptides have in the last years clarified several issues concerning peptide folding equilibria and interpretation of experimental data, especially from NMR and CD spectroscopy. These simulations involved different temperatures, pH-values, ionic strengths, solvents, and force-field parameters, but a variation of these factors for one beta-peptide has not yet been done. To investigate the influence of varying these factors, we analyze the helix stability of an all-beta3-icosapeptide bearing all 20 proteinogenic amino acid side chains, which is experimentally observed to fold into a 3(14)-helix in methanol but not in water. Structural aspects, such as hydrogen-bonded rings and salt bridges, are discussed and a comparison with NMR primary (NOE distance bounds and 3J-values) and secondary (NMR derived model structures) data is made. We further investigate the reasons for the 3(14)-helix stability/instability in methanol/water. Of all factors studied, the presence of counterions seems to be the one inducing most significant effects in the simulations., (Copyright 2006 Wiley Periodicals, Inc.)

Molecular dynamic simulations have been performed for wild-type Hydrogenobacter thermophilus cytochrome c(552), a b-type variant of the protein, and the apo state with the heme prosthetic group removed. In the b-type variant, Cys 10 and Cys 13 were mutated to alanine residues, and so the heme group was no longer covalently bound to the protein. Two 8-ns simulations have been performed for each system at 298 and 360 K. The simulations of the wild-type protein at 298 K show a very close agreement with experimental NMR data. A fluxional process involving the side chain of Met 59, which coordinates to the heme iron, is observed in accord with proposals from NMR studies. Overall, the structure and dynamical behavior of the protein during the simulations of the b-type variant is closely similar to that of the wild-type protein. However, side chains in the heme-binding site show larger fluctuations in the b-type variant simulation at 360 K. In addition, structural changes are seen for a number of residues close to the heme group, particularly Gly 22 and Ser 51. The simulations of the apo state show significant conformational changes for residues 50-59. These residues form a loop region, which packs over the heme group in the wild-type protein and hydrogen bonds to the heme propionate groups. In the absence of heme, in the apo state simulations, these residues form short but persistent regions of beta-sheet secondary structure. These could provide nucleation sites for the conversion to amyloid fibrils., ((c) 2006 Wiley-Liss, Inc.)

Cyclophilins are proteins that catalyze X-proline cis-trans interconversion, where X represents any amino acid. Its mechanism of action has been investigated over the past years but still generates discussion, especially because until recently structures of the ligand in the cis and trans conformations for the same system were lacking. X-ray crystallographic structures for the complex cyclophilin A and HIV-1 capsid mutants with ligands in the cis and trans conformations suggest a mechanism where the N-terminal portion of the ligand rotates during the cis-trans isomerization. However, a few years before, a C-terminal rotating ligand was proposed to explain NMR solution data. In the present study we use molecular dynamics (MD) simulations to generate a trans structure starting from the cis structure. From simulations starting from the cis and trans structures obtained through the rotational pathways, the seeming contradiction between the two sets of experimental data could be resolved. The simulated N-terminal rotated trans structure shows good agreement with the equivalent crystal structure and, moreover, is consistent with the NMR data. These results illustrate the use of MD simulation at atomic resolution to model structural transitions and to interpret experimental data.

Two designed ankyrin repeat (AR) proteins (E3_5 and E3_19) are high homologous (with about 87% sequence identity) and their crystal structures have a Calpha atom-positional root-mean-square difference of about 0.14 nm. However, it was found that E3_5 is considerably more stable than E3_19 in guanidinium hydrochloride and thermal denaturation experiments. With the goal of providing insights into the various factors contributing to the stabilities of the designed AR proteins and suggesting possible mutations to enhance their stabilities, homology modeling and molecular dynamics (MD) simulations with explicit solvent have been performed. Because the crystal structure of E3_19 was solved later than that of E3_5, a homology model of E3_19 based on the crystal structure of E3_5 was also used in the simulations. E3_5 shows a very stable trajectory in both crystal and solution simulations. In contrast, the C-terminal repeat of E3_19 unfolds in the simulations starting from either the modeled structure or the crystal structure, although it has a sequence identical to that of E3_5. A continuum electrostatic model was used to estimate the effect of single mutations on protein stability and to study the interaction between the internal ARs and the C-terminal capping AR. Mutations involving charged residues were found to have large effects on stability. Due to the difference in charge distribution in the internal ARs of E3_19 and E3_5, their interaction with the C-terminal capping AR is less favorable in E3_19. The simulation trajectories suggest that the stability of the designed AR proteins can be increased by optimizing the electrostatic interactions within and between the different repeats., ((c) 2006 Wiley-Liss, Inc.)

Experimental nuclear magnetic resonance results for the Arc Repressor have shown that this dimeric protein dissociates into a molten globule at high pressure. This structural change is accompanied by a modification of the hydrogen-bonding pattern of the intermolecular beta-sheet: it changes its character from intermolecular to intramolecular with respect to the two monomers. Molecular dynamics simulations of the Arc Repressor, as a monomer and a dimer, at elevated pressure have been performed with the aim to study this hypothesis and to identify the major structural and dynamical changes of the protein under such conditions. The monomer appears less stable than the dimer. However, the complete dissociation has not been seen because of the long timescale needed to observe this phenomenon. In fact, the protein structure altered very little when increasing the pressure. It became slightly compressed and the dynamics of the side-chains and the unfolding process slowed down. Increasing both, temperature and pressure, a tendency of conversion of intermolecular into intramolecular hydrogen bonds in the beta-sheet region has been detected, supporting the mentioned hypothesis. Also, the onset of denaturation of the separated chains was observed., (Proteins 2006. (c) 2006 Wiley-Liss, Inc.)

A polarizable model for simulation of liquid methanol, compatible with the COS/G2 water model, has been developed using the Charge-on-Spring (COS) technique. The model consists of three point charges, with one polarizable center on the oxygen atom. The Lennard-Jones parameters on the oxygen atom together with the molecular polarizability were varied to reproduce the experimental heat of vaporization and density of liquid methanol at ambient conditions. We examined the energies of various methanol dimers in the gas phase and compared them with values obtained from ab initio calculations. The model was then used to study the thermodynamic, dynamic, structural, and dielectric properties of liquid methanol as well as of a methanol-water mixture. A microscopic picture of the structure of pure liquid methanol and of the methanol-water mixture is provided. Good agreement was found between the results from our model simulations and available experimental and ab initio calculation data. In particular, the experimental dielectric permittivity of 32 could be reproduced, which had been shown to be difficult when using nonpolarizable models., (Copyright 2006 Wiley Periodicals, Inc.)

Binding of a small molecule to a macromolecular target reduces its conformational freedom, resulting in a negative entropy change that opposes the binding. The goal of this study is to estimate the configurational entropy change of two minor-groove-binding ligands, netropsin and distamycin, upon binding to the DNA duplex d(CGCGAAAAACGCG).d(CGCGTTTTTCGCG). Configurational entropy upper bounds based on 10-ns molecular dynamics simulations of netropsin and distamycin in solution and in complex with DNA in solution were estimated using the covariance matrix of atom-positional fluctuations. The results suggest that netropsin and distamycin lose a significant amount of configurational entropy upon binding to the DNA minor groove. The estimated changes in configurational entropy for netropsin and distamycin are -127 J K(-1) mol(-1) and -104 J K(-1) mol(-1), respectively. Estimates of the configurational entropy contributions of parts of the ligands are presented, showing that the loss of configurational entropy is comparatively more pronounced for the flexible tails than for the relatively rigid central body.

Single-chain and single-fragment configurational entropies of lipid tails in hydrated lipid bilayers are evaluated from molecular dynamics simulations using the quasi-harmonic approximation. The entropy distribution along individual acyl tails is obtained and compared to that of corresponding hydrocarbon chains in the liquid phase. We consider pure dipalmitoylphosphatidylcholine and mixed dioleoylphosphatidylcholine/dioleoylphosphatidylethanolamine bilayers. The systems are modeled at different levels of spatial resolution: In an atomic-level (AL) model all (heavy) atoms are explicitly simulated; in a coarse-grained (CG) model particles (beads) representing groups of covalently bound atoms are used, which map approximately four non-hydrogen atoms to one interaction site. Single-chain and single-fragment entropies and correlations between the motions of (single) acyl chains are compared. A good correspondence is found between the flexibility of the AL and CG models. The loss in configurational entropy due to the reduction in the number of degrees of freedom upon coarse-graining of the model is estimated. The CG model shows about 4 times faster convergence of the chain entropies than the more detailed AL model. Corrections to the quasi-harmonic entropy estimates were found to be small for the CG model. For the AL model, the correction due to mode anharmonicities is small, but the correction due to pairwise (supralinear) mode correlations is sizable.

Computation based on molecular models is playing an increasingly important role in biology, biological chemistry, and biophysics. Since only a very limited number of properties of biomolecular systems is actually accessible to measurement by experimental means, computer simulation can complement experiment by providing not only averages, but also distributions and time series of any definable quantity, for example, conformational distributions or interactions between parts of systems. Present day biomolecular modeling is limited in its application by four main problems: 1) the force-field problem, 2) the search (sampling) problem, 3) the ensemble (sampling) problem, and 4) the experimental problem. These four problems are discussed and illustrated by practical examples. Perspectives are also outlined for pushing forward the limitations of biomolecular modeling.

Explicit-solvent molecular dynamics simulations (50 ns, 300 K) of the eight reducing glucose disaccharides (kojibiose, sophorose, nigerose, laminarabiose, maltose, cellobiose, isomaltose, and gentiobiose) have been carried out using the GROMOS 45A4 force field (including a recently reoptimized carbohydrate parameter set), to investigate and compare their conformational preferences, intramolecular hydrogen-bonding patterns, torsional dynamics, and configurational entropies. The calculated average values of the glycosidic torsional angles agree well with available experimental data, providing validation for the force field and simulation methodology employed in this study. These simulations show in particular that: 1) (1-->6)-linked disaccharides are characterized by an increased flexibility, the absence of any persistent intramolecular hydrogen bond and a significantly higher configurational entropy (compared to the other disaccharides); 2) cellobiose presents a highly persistent interresidue hydrogen bond and a significantly lower configurational entropy (compared to the other disaccharides); 3) persistent hydrogen bonds are observed for all disaccharides (except (1-->6)-linked) and typically involve a hydrogen donor in the reducing residue and an acceptor in the nonreducing one; 4) the probability distributions associated with the glycosidic dihedral angles and psi are essentially unimodal for all disaccharides, and full rotation around these angles occurs at most once or twice for (never for psi) on the 50-ns timescale; and 5) the timescales associated with torsional transitions (except around and psi) range from approximately 30 ps (rotation of hydroxyl groups) to the nanosecond range (rotation of the lactol and hydroxymethyl groups, and around the omega-glycosidic dihedral angle in (1-->6)-linked disaccharides).

The relative stabilities of homochiral and heterochiral forms of selected dipeptides, AA, AS, AC, AV, AF, AD, AK, tripeptides, AAA, AVA, and an acetylpentapeptide, AcGLSFA, have been calculated using thermodynamic integration protocols and the GROMOS 53A6 force field. Integration pathways have been designed that produce minimal disturbance to the system, including the use of soft atoms, low-energy intermediates, and chiral inversion of the smaller amino acid in the peptide. Comparison of the results obtained by thermodynamic integration between the diastereomeric forms (in explicit water, at 300 K) and from exhaustive global minimum-energy searches for the individual dipeptides (implicit water, epsilon = 78, 0 K) suggests that entropic contributions to the relative stability of the chiral forms are important. This conclusion is supported by the results of explicit calculation of the effect of temperature on the relative stability of alanylvalylalanine diastereomers. The Gibbs free energy calculations predict that at ambient temperature and pressure homochiral dipeptides with small side chains or polar groups in the vicinity of the peptide backbone, AA, AS, and AD, are more stable than their heterochiral counterparts by fractions of a kJ/mol. For bigger side chains, AC, AV, AF, and AK, the heterochiral diastereomers appear to be more stable. Predicted relative stabilities are in line with observations reported in the literature for AE and YY. Excellent agreement is found for the calculated and experimentally determined relative stabilities of the diastereomers of the dipeptide AA and of all-L AcGLSFA and its diastereomer containing D-serine in the central position. Addition of counterions to the solvent box has no significant effects on charged and neutral forms. From the present findings it would appear unlikely that the intrinsic stability difference between homo- and heterochiral dipeptides has been a driving force in a primordial selection process leading to the incorporation of amino acids with a single enantiomeric configuration in natural proteins., ((c) 2006 Wiley Periodicals, Inc.)

A method to enhance sampling of rare events is presented. It makes use of distance or dihedral-angle restraints to overcome an energy barrier separating two metastable states or to stabilize a transition state between the two metastable states. In order not to perturb these metastable end states themselves, a prefactor is introduced into the restraining energy function, which smoothly increases the weight of this function from zero to one at the transition state or on top of the separating energy barrier and then decreases the weight again to zero at the final state. The method is combined with multi-configurational thermodynamic integration and applied to two biomolecular systems, which were difficult to treat using standard thermodynamic integration. As first example the free energy difference of a cyclic alpha-aminoxy-hexapeptide-ion complex upon changing the ion from Cl- to Na+ was calculated. A large conformational rearrangement of the peptide was necessary to accommodate this change. Stabilizing the transition state by (hidden) restraints facilitates that. As a second example, the free energy difference between the 4C1 and the 1C4 conformation of beta-D-glucopyranoside was calculated. In unrestrained simulations the change from the 4C1 into the 1C4 conformation was never observed because of the high energy barrier separating the two states. Using (hidden) restraints, the transition from the 4C1 into the 1C4 state and back could be enforced without perturbing the end states. As comparison, for the same transitions the potential of mean force as obtained by using dihedral-angle constraints is provided.

Molecular liquids can be modeled at different levels of spatial resolution. In atomic-level (AL) models, all (heavy) atoms can be explicitly simulated. In coarse-grained (CG) models, particles (beads) that represent groups of covalently bound atoms are used as elementary units. Ideally, a CG model should reproduce the thermodynamic and structural properties of the corresponding AL model after mapping to the lower-resolution scale. In the present work, two such models are investigated: (i) the classical GROMOS atomic-level model; (ii) a CG model recently proposed by Marrink et al., which maps approximately four non-hydrogen atoms to one bead [J. Phys. Chem. B 2004, 108, 750]. The study is restricted to n-alkanes whose aliphatic fragments are abundantly found in lipids of biological interest. Additionally, cis-9-octadecene is included, as a template chain of the lipid dioleoylphosphatidylcholine (DOPC). The two representations of molecules in the liquid phase are compared in terms of average molecular structures, extent of configurational space sampled, and single-molecule entropies. An approximate method is used to estimate the rotational contributions to the absolute configurational entropy. Good correspondence between the AL and CG representations is found. The loss in configurational entropy due to the reduction in degrees of freedom upon coarse-graining of the model is estimated.

A method to combine fine-grained and coarse-grained simulations is presented. The coarse-grained particles are described as virtual particles defined by the underlying fine-grained particles are described as virtual particles defined by the underlying fine-grained particles. The contribution of the two grain levels to the interaction between particles is specified by a grain-level parameter lambda. Setting lambda = 0 results in a completely fine-grained simulation, whereas lambda = 1 yields a simulation governed by the coarse-grained potential energy surface with small contributions to keep the fine-grained covalently bound particles together. Simulations at different lambda values may be coupled using the replica-exchange molecular dynamics method to achieve enhanced sampling at the fine-grained level.

It has been suggested that the stability of a beta-peptide helical fold is affected by the interplay between the electrical charge of terminal groups and the dipole due to the helical conformation, the so-called charge-dipole stabilization; the numerical simulations presented herein test that suggestion. The motions of two beta-peptide oligomers, each of which has been shown by NMR spectroscopy to fold into a different helical conformation, have been simulated. The simulated motions bear out empirical observations as to the effect of chemical protection of terminal groups on the stability of beta-peptide helical folds and they support the hypothesis of charge-dipole stabilization., (2006 Wiley-Liss, Inc.)

Simulated molecular dynamics trajectories of proteins and nucleic acids are often compared with nuclear magnetic resonance (NMR) data for the purposes of assessing the quality of the force field used or, equally important, trying to interpret ambiguous experimental data. In particular, nuclear Overhauser enhancement (NOE) intensities or atom-atom distances derived from them are frequently calculated from the simulated ensembles because the distance restraints derived from NOEs are the key ingredient in NMR-based protein structure determination. In this study, we ask how diverse and nonnative-like an ensemble of structures can be and still match the experimental NOE distance upper bounds well. We present two examples in which simulated ensembles of highly nonnative polypeptide structures (an unfolded state ensemble of the villin headpiece and a high-temperature denatured ensemble of lysozyme) are shown to match fairly well the experimental NOE distance upper bounds from which the corresponding native structures were derived. For example, the unfolded ensemble of villin headpiece, which is on average 0.90 +/- 0.13 nm root-mean-square deviation away from the native NMR structure, deviates from the experimental restraints by only 0.027 nm on average. However, this artificially good agreement is largely a consequence of 1) the highly nonlinear effects of r(-6) (or r(-3)) averaging and 2) focusing only on the experimentally observed set of NOE bounds. Namely, in addition to the experimentally observed NOEs, both simulated ensembles (especially the villin ensemble) also predict a large number of NOEs, which are not seen in the experiment. If these are taken into account, the agreement between simulation and experiment gets markedly worse, as it should, given the nonnative nature of the underlying simulated ensembles. In light of the examples given, we conclude that comparing experimental NOE distance restraints with large simulated ensembles provides just by itself only limited information about the quality of simulation., (2006 Wiley-Liss, Inc.)

Recently, the GROMOS biomolecular force field parameter set 53A6--which has been parametrized to reproduce experimentally determined free enthalpies of hydration and solvation in cyclohexane of amino acid side-chain analogs--was presented. To investigate the transferability of the new parameter set, we calculated free enthalpies of solvation of a range of polar and apolar compounds in different solvents (methanol, dimethyl sulfoxide (DMSO), acetonitrile, and acetone) from molecular dynamics simulations using the GROMOS 53A6 force field. For methanol and DMSO, parameters were used that are available in the 53A6 parameter set. For acetonitrile, a recently developed model was taken and for acetone, two models available in literature were used. We found that trends in and values for the solvation free enthalpies are in satisfactory agreement with experiment, except for the solvation in acetone for which deviations from experiment can be explained in terms of the properties of the models used.

The structural properties of melittin, a small amphipathic peptide found in the bee venom, are investigated in three different environments by molecular dynamics simulation. Long simulations have been performed for monomeric melittin solvated in water, in methanol, and shorter ones for melittin inserted in a dimyristoylphosphatidylcholine bilayer. The resulting trajectories were analysed in terms of structural properties of the peptide and compared to the available NMR data. While in water and methanol solution melittin is observed to partly unfold, the peptide retains its structure when embedded in a lipid bilayer. The latter simulation shows good agreement with the experimentally derived (3)J-coupling constants. Generally, it appears that higher the stability of the helical conformation of melittin, lower is the dielectric permittivity of the environment. In addition, peptide-lipid interactions were investigated showing that the C-terminus of the peptide provides an anchor to the lipid bilayer by forming hydrogen bonds with the lipid head groups.

The morphogenic opiate pentapeptide leucine-enkephalin (lenk) in a hydrated dimyristoylphosphatidylcholine (DMPC) bilayer is studied using NMR spectroscopy and molecular dynamics simulation. Contrary to the frequent assumption that the peptide attains a single fixed conformation in the presence of membranes, we find that the lenk molecule is flexible, switching between specific bent conformations. The constraints to the orientation of the aromatic rings that are identified by the NMR experiment are found by the MD simulation to be related to the depth of the peptide in the bilayer. The motion of the N-H vectors of the peptide bonds with respect to the magnetic field direction as observed by MD largely explain the magnitude of the observed residual dipolar coupling (RDC), which are much reduced over the static (15)N-(1)H coupling. The measured RDCs are nevertheless significantly larger than the predicted ones, possibly due the absence of long-time motions in the simulations. The conformational behavior of lenk at the DMPC surface is compared to that in the aqueous solution, both in the neutral and in the zwitterionic forms.

This report examines the effect of a decrease in solvent viscosity on the simulated folding behaviour of a beta-peptide heptamer in methanol. Simulations of the molecular dynamics of the heptamer H-beta3-HVal-beta3-HAla-beta3-HLeu-(S,S)-beta3-HAla(alphaMe)-beta3-HVal-beta3-HAla-beta3-HLeu-OH in methanol, with an explicit representation of the methanol molecules, were performed for 80 ns at various solvent viscosities. The simulations indicate that at a solvent viscosity of one third of that of methanol, only the dynamic aspects of the folding process are altered, and that the rate of folding is increased. At a viscosity of one tenth of that of methanol, insufficient statistics are obtained within the 80 ns period. We suggest that 80 ns is an insufficient time to reach conformational equilibrium at very low viscosity because the dependence of the folding rate of a beta-peptide on solvent viscosity has two regimes; a result that was observed in another computational study for alpha-peptides.

The influence of charged side chains on the folding-unfolding equilibrium of beta-peptides was investigated by means of molecular dynamics simulations. Four different peptides containing only negatively charged side chains, positively charged side chains, both types of charged side chains (with the ability to form stabilizing salt bridges) or no charged side chains were studied under various conditions (different simulation temperatures, starting structures and solvent environment). The NMR solution structure in methanol of one of the peptides (A) has already been published; the synthesis and NMR analysis of another peptide (B) is described here. The other peptides (C and D) studied herein have hitherto not been synthesized. All four peptides A-D are expected to adopt a left-handed 3(14)-helix in solution as well as in the simulations. The resulting ensembles of structures were analyzed in terms of conformational space sampled by the peptides, folding behavior, structural properties such as hydrogen bonding, side chain-side chain and side chain-backbone interactions and in terms of the level of agreement with the NMR data available for two of the peptides. It was found that the presence of charged side chains significantly slows down the folding process in methanol solution due to the stabilization of intermediate conformers with side chain-backbone interactions. In water, where the solvent competes with the solute-solute polar interactions, the folding process to the 3(14)-helix is faster in the simulations.

We present the latest version of the Groningen Molecular Simulation program package, GROMOS05. It has been developed for the dynamical modelling of (bio)molecules using the methods of molecular dynamics, stochastic dynamics, and energy minimization. An overview of GROMOS05 is given, highlighting features not present in the last major release, GROMOS96. The organization of the program package is outlined and the included analysis package GROMOS++ is described. Finally, some applications illustrating the various available functionalities are presented., ((c) 2005 Wiley Periodicals, Inc.)

NMR is one of the most used techniques to resolve structure of proteins and peptides in solution. However, inconsistencies may occur due to the fact that a polypeptide may adopt more than one conformation. Since the NOE distance bounds and (3)J-values used in such structure determination represent a nonlinear average over the total ensemble of conformers, imposition of NOE or (3)J-value restraints to obtain one unique conformation is not an appropriate procedure in such cases. Here, we show that unrestrained MD simulation of a solute in solution using a high-quality force field yields a conformational ensemble that is largely compatible with the experimental NMR data on the solute. Four 100 ns MD simulations of two forms of a nine-residue beta-peptide in methanol at two temperatures produced conformational ensembles that were used to interpret the NMR data on this molecule and resolve inconsistencies between the experimental NOEs. The protected and unprotected forms of the beta-peptide adopt predominantly a 12/10-helix in agreement with the qualitative interpretation of the NMR data. However, a particular NOE was not compatible with this helix indicating the presence of other conformations. The simulations showed that 3(14)()-helical structures were present in the ensemble of the unprotected form and that their presence correlates with the fulfillment of the particular NOE. Additionally, all inter-hydrogen distances were calculated to compare NOEs predicted by the simulations to the ones observed experimentally. The MD conformational ensembles allowed for a detailed and consistent interpretation of the experimental data and showed the small but specific conformational differences between the protected and unprotected forms of the peptide.

The free energies of hydration for ammonia and mono-, di-, and trimethylated amines experimentally show an unexpected trend that has, in the past, been difficult to reproduce computationally. Absolute and relative free energies of hydration of these compounds were calculated using the OPLS all-atom and the united-atom GROMOS force fields. Both force fields reproduce the relative free energy of hydration, but the absolute free energies of hydration were only reproduced within kBT when using the recently developed GROMOS parameter set 53A6. Relative free energies of solvation in chloroform could also be reproduced indicating a proper partitioning of the compounds between polar and apolar media. Overall we conclude that it is possible to resolve the amine hydration problem using a simple united-atom force field.

Polyproline type II (PPII) helix has emerged recently as the dominant paradigm for describing the conformation of unfolded polypeptides. However, most experimental observables used to characterize unfolded proteins typically provide only short-range, sequence-local structural information that is both time- and ensemble-averaged, giving limited detail about the long-range structure of the chain. Here, we report a study of a long-range property: the radius of gyration of an alanine-based peptide, Ace-(diaminobutyric acid)2-(Ala)7-(ornithine)2-NH2. This molecule has previously been studied as a model for the unfolded state of proteins under folding conditions and is believed to adopt a PPII fold based on short-range techniques such as NMR and CD. By using synchrotron radiation and small-angle x-ray scattering, we have determined the radius of gyration of this peptide to be 7.4 +/- 0.5 angstroms, which is significantly less than the value expected from an ideal PPII helix in solution (13.1 angstroms). To further study this contradiction, we have used molecular dynamics simulations using six variants of the AMBER force field and the GROMOS 53A6 force field. However, in all cases, the simulated ensembles underestimate the PPII content while overestimating the experimental radius of gyration. The conformational model that we propose, based on our small angle x-ray scattering results and what is known about this molecule from before, is that of a very flexible, fluctuating structure that on the level of individual residues explores a wide basin around the ideal PPII geometry but is never, or only rarely, in the ideal extended PPII helical conformation.

We report the investigation of two 16-residue peptides in aqueous solution by means of molecular-dynamics simulations. The peptides constitute the C- and N-terminal halves of the 33-residue monomer whose dimer constitutes the leucine zipper of the yeast transcriptional activator, denoted GCN4-p1. To examine a hypothesis about coiled-coil formation, in which the C-terminal half contains a helix-formation trigger site absent in the N-terminal half, experimental studies of the two peptides have determined their helix propensities under several conditions of temperature, pH, and salt concentration with circular dichroism. An NMR experiment provides additional evidence. At temperatures of 278 and 325 K and pH 7.5, mixtures of alpha- and pi-helical secondary structure constitute the most probable conformations in both C- and N-terminal halves. A bifurcated salt bridge between Arg25 and Glu22/20 correlates with the structural fluctuations of the C-terminal half. It also exhibits a persistent loop at the N-terminal end involving the side chains of His18 and Glu22, which is reminiscent of helix-capping boxes. Nonreversible unfolding appears to occur abruptly in the Arg25 mutant, suggesting a cooperative event. Analysis does not indicate that the N-terminal half is less stable than the C-terminal half, indicating that 100 ns is too short a period to observe complete unfolding.

Through the use of the one-step perturbation approach, 130 free energies of base stacking and 1024 free energies of base pairing in DNA have been calculated from only five simulations of a nonphysical reference state. From analysis of a diverse set of 23 natural and unnatural bases, it appears that stacking free energies and stacking conformations play an important role in pairing of DNA nucleotides. On the one hand, favourable pairing free energies were found for bases that do not have the possibility to form canonical hydrogen bonds, while on the other hand, good hydrogen-bonding possibilities do not guarantee a favourable pairing free energy if the stacking of the bases dictates an unfavourable conformation. In this application, the one-step perturbation approach yields a wealth of both energetic and structural information at minimal computational cost.