Search

In protein simulation or structure refinement based on values of observable quantities measured in (aqueous) solution, solvent (water) molecules may be explicitly treated, omitted, or represented by a potential of mean-solvation-force term, depending on protein coordinates only, in the force field used. These three approaches are compared for hen egg white lysozyme (HEWL). This 129-residue non-spherical protein contains a variety of secondary-structure elements, and ample experimental data are available: 1630 atom-atom Nuclear Overhauser Enhancement (NOE) upper distance bounds, 213 3  J-couplings and 200 S 2 order parameters. These data are used to compare the performance of the three approaches. It is found that a molecular dynamics (MD) simulation in explicit water approximates the experimental data much better than stochastic dynamics (SD) simulation in vacuo without or with a solvent-accessible-surface-area (SASA) implicit-solvation term added to the force field. This is due to the missing energetic and entropic contributions and hydrogen-bonding capacities of the water molecules and the missing dielectric screening effect of this high-permittivity solvent. Omission of explicit water molecules leads to compaction of the protein, an increased internal strain, distortion of exposed loop and turn regions and excessive intra-protein hydrogen bonding. As a consequence, the conformation and dynamics of groups on the surface of the protein, which may play a key role in protein-protein interactions or ligand or substrate binding, may be incorrectly modelled. It is thus recommended to include water molecules explicitly in structure refinement of proteins in aqueous solution based on nuclear magnetic resonance (NMR) or other experimentally measured data., (© 2022. The Author(s).)

The upper and lower lip regions in lysozyme from bacteriophage lambda (λ-lysozyme) are flexible in solution and exhibit two different conformations in crystal structures of the protein. MD simulations have been used to characterize the structure and dynamics of these lip regions, which surround the active site. Ten different simulations have been run including those with restraining to experimental NOE distance and (1)H-(15)N order parameter data. The simulations show that the lower lip region, although undergoing considerable backbone fluctuations, contains two persistent β-strands. In the upper lip region, a wide range of conformations are populated and it is not clear from the available data whether some helical secondary structure is present. The work provides a clear example of the advantages of combining MD simulations with experimental data to obtain a structural interpretation of the latter. In this case, time-averaged order parameter restraining has played an essential role in enabling convergence between two different starting structures and identifying the extent to which flexible regions in solution can contain persistent secondary structure.

The method of one-step perturbation can be used to predict from a single molecular dynamics simulation the values of observable quantities as functions of variations in the parameters of the Hamiltonian or biomolecular force field used in the simulation. The method is used to predict violations of nuclear overhauser effect (NOE) distance bounds measured in nuclear magnetic resonance (NMR) experiments by atom-atom distances of the NOE atom pairs when varying force-field parameters. Predictions of NOE distance bound violations between different versions of the GROMOS force field for a hexa-β-peptide in solution show that the technique works for rather large force-field parameter changes as well as for very different NOE bound violation patterns. The effect of changing individual force-field parameters on the NOE distance bound violations of the β-peptide and an α-peptide was investigated too. One-step perturbation, which in this case is equivalent to reweighting configurations, constitutes an efficient technique to predict many values of different quantities from a single conformational ensemble for a particular system, which makes it a powerful force-field development technique that easily reduces the number of required separate simulations by an order of magnitude.

It is often assumed that intramolecular hydrogen-bonding (H-bonding) exerts a significant influence on the conformational properties of aqueous (bio-)polymers. To discuss this statement, one should, however, distinguish between solvent-exposed and buried H-bonds, and between their respective roles in promoting stability (i.e., as a driving force) and specificity (for which the term steering force is introduced here). In this study, the role of solvent-exposed H-bonding in carbohydrates as a driving or steering force is probed using explicit-solvent molecular dynamics simulations with local elevation umbrella sampling in the simple context of cellobiose stereoisomers. More specifically, four β(1→4)-linked D-aldohexopyranose disaccharides are considered, which present a different stereochemisty of the potentially H-bonding groups neighboring the glycosidic linkage. Although the epimerization may largely alter the intramolecular trans-glycosidic H-bonding pattern, it is found to have only very limited influence on the Ramachandran free-energy map of the disaccharide, a loss of intramolecular H-bonding being merely compensated for by an enhancement of the interaction with the solvent molecules. This finding suggests that solvent-exposed trans-glycosidic H-bonding (and in particular the HO'(3)→O5 H-bond) is not the cause of the 21-helical secondary structure characteristic of cellooligosaccharides, but rather the opportunistic consequence of a sterically and stereoelectronically dictated conformational preference. In other words, for these compounds, solvent-exposed H-bonding appears to represent a minor (possibly adverse) conformational driving as well as steering force.

Simulation of the dynamics of a protein in aqueous solution using an atomic model for both the protein and the many water molecules is still computationally extremely demanding considering the time scale of protein motions. The use of supra-atomic or supra-molecular coarse-grained (CG) models may enhance the computational efficiency, but inevitably at the cost of reduced accuracy. Coarse-graining solvent degrees of freedom is likely to yield a favourable balance between reduced accuracy and enhanced computational speed. Here, the use of a supra-molecular coarse-grained water model that largely preserves the thermodynamic and dielectric properties of atomic level fine-grained (FG) water in molecular dynamics simulations of an atomic model for four proteins is investigated. The results of using an FG, a CG, an implicit, or a vacuum solvent environment of the four proteins are compared, and for hen egg-white lysozyme a comparison to NMR data is made. The mixed-grained simulations do not show large differences compared to the FG atomic level simulations, apart from an increased tendency to form hydrogen bonds between long side chains, which is due to the reduced ability of the supra-molecular CG beads that represent five FG water molecules to make solvent-protein hydrogen bonds. But, the mixed-grained simulations are at least an order of magnitude faster than the atomic level ones.

Deriving structural information about a protein from NMR experimental data is still a non-trivial challenge to computational biochemistry. This is because of the low ratio of the number of independent observables to the number of molecular degrees of freedom, the approximations involved in the different relationships between particular observable quantities and molecular conformation, and the averaged character of the experimental data. For example, protein (3)J-coupling data are seldom used for structure refinement because of the multiple-valuedness and limited accuracy of the Karplus relationship linking a (3)J-coupling to a torsional angle. Moreover, sampling of the large conformational space is still problematic. Using the 99-residue protein plastocyanin as an example we investigated whether use of a thermodynamically calibrated force field, inclusion of solvent degrees of freedom, and application of adaptive local-elevation sampling that accounts for conformational averaging produces a more realistic representation of the ensemble of protein conformations than standard single-structure refinement in a non-explicit solvent using restraints that do not account for averaging and are partly based on non-observed data. Yielding better agreement with observed experimental data, the protein conformational ensemble is less restricted than when using standard single-structure refinement techniques, which are likely to yield a picture of the protein which is too rigid.

New parameter sets of the GROMOS biomolecular force field, 54A7 and 54B7, are introduced. These parameter sets summarise some previously published force field modifications: The 53A6 helical propensities are corrected through new φ/ψ torsional angle terms and a modification of the N-H, C=O repulsion, a new atom type for a charged -CH(3) in the choline moiety is added, the Na(+) and Cl(-) ions are modified to reproduce the free energy of hydration, and additional improper torsional angle types for free energy calculations involving a chirality change are introduced. The new helical propensity modification is tested using the benchmark proteins hen egg-white lysozyme, fox1 RNA binding domain, chorismate mutase and the GCN4-p1 peptide. The stability of the proteins is improved in comparison with the 53A6 force field, and good agreement with a range of primary experimental data is obtained.

The bacterial outer membrane protein OmpX from Escherichia coli has been investigated by molecular dynamics simulations when embedded in a phospholipid bilayer and as a protein-micelle aggregate. The resulting simulation trajectories were analysed in terms of structural and dynamic properties of the membrane protein. In agreement with experimental observations, highest relative stability was found for the β-barrel region that is embedded in the lipophilic phase, whereas an extracellular protruding β-sheet, which is a unique structural feature of OmpX that supposedly plays an important role in cell adhesion and invasion, shows larger structure fluctuations. Additionally, we investigated water permeation into the core of the β-barrel protein, which contains a tight salt-bridge and hydrogen-bond network, so that extensive water flux is unlikely. Differences between the bilayer and the micellar system were observed in the length of the barrel and its position inside the lipid environment, and in the protein interactions with the hydrophilic part of the lipids near the lipid/water interface. Those variations suggest that micelles and other detergent environments might not offer a wholly membrane-like milieu to promote adoption of the physiological conformational state by OmpX.

We have studied two different beta-peptides in methanol using explicit solvent molecular dynamics simulations and the GROMOS 53A6 force field: a heptapeptide (peptide 1) expected to form a left-handed 3(14)-helix, and a hexapeptide (peptide 2) expected to form a beta-hairpin in solution. Our analysis has focused on identifying and analyzing the stability of the dominant secondary structure conformations adopted by the peptides, as well as on comparing the experimental NOE distance upper bounds and 3J-coupling values with their counterparts calculated on the basis of the simulated ensembles. Moreover, we have critically compared the present results with the analogous results obtained with the GROMOS 45A3 (peptide 1) and 43A1 (peptide 2) force fields. We conclude that within the limits of conformational sampling employed here, the GROMOS 53A6 force field satisfactorily reproduces experimental findings regarding the behavior of short beta-peptides, with accuracy that is comparable to but not exceeding that of the previous versions of the force field. GCE legend Conformational clustering analysis of the simulated ensemble of a ss-hexapeptide with two different simulation setups (a and b). The central members of all of the clusters populating more than 5% of all of the structures are shown, together with the most dominant hydrogen bonds and the corresponding percentages of cluster members containing them.

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 quality of biomolecular dynamics simulations relies critically on the force field that is used to describe the interactions between particles in the system. Force fields, which are generally parameterized using experimental data on small molecules, can only prove themselves in realistic simulations of relevant biomolecular systems. In this work, we begin the validation of the new 53A6 GROMOS parameter set by examining three test cases. Simulations of the well-studied 129 residue protein hen egg-white lysozyme, of the DNA dodecamer d(CGCGAATTCGCG)(2), and a proteinogenic beta(3)-dodecapeptide were performed and analysed. It was found that the new parameter set performs as well as the previous parameter sets in terms of protein (45A3) and DNA (45A4) stability and that it is better at describing the folding-unfolding balance of the peptide. The latter is a property that is directly associated with the free enthalpy of hydration, to which the 53A6 parameter set was parameterized.

Circular dichroism spectra of two beta-peptides, i.e. peptides composed of beta-amino acids, calculated using ensembles of configurations obtained by molecular dynamics simulation are presented. The calculations were based on 200 ns simulations of a beta-heptapeptide in methanol at 298 K and 340 K and a 50 ns simulation of a beta-hexapeptide in methanol at 340 K. In the simulations the peptides sampled both folded (helical) and unfolded structures. Trajectory structures with common backbone conformations were identified and grouped into clusters. The CD spectra were calculated for individual structures, based on peptide-group dipole transition moments obtained from semi-empirical molecular orbital theory and using the so-called matrix method. The single-structure spectra were then averaged over entire trajectories and over clusters of structures. Although certain features of the experimental CD spectra of the beta-peptides are reproduced by the trajectory-average spectra, there exist clear differences between the two sets of spectra in both wavelength and peak intensities. The analysis of individual contributions to the average spectra shows that, in general, the interpretation of a CD signal in terms of a single structure is not possible. Moreover, there is a large variation in the CD spectra calculated for a set of individual structures that belong to the same cluster, even when a structurally tight clustering criterion is used. This indicates that the CD spectra of these peptides are very sensitive to small local structural differences.

The performance of the GROMOS96 parameter set 45A3 developed for aliphatic alkanes is tested on a bilayer of dipalmitoylphosphatidylcholine (DPPC) in water in the liquid-crystalline L(alpha) phase. Variants of the force-field parameter set as well as different sets of simulation conditions or simulation parameter sets are evaluated. In the case of the force-field parameters, the van der Waals constants for the non-bonded interaction of the ester carbonyl carbon and the partial charges and charge group definition of the phosphatidylcholine head group are examined. On the methodological side, different cut-off distances for the non-bonded interactions, use of a reaction-field force due to long-range electrostatic interactions, the frequency of removal of the centre of mass motion and the strength of the coupling of the pressure of the system to the pressure bath are tested. The area per lipid, as a measure of structure, the order parameters of the chain carbons, as a measure of membrane fluidity, and the translational diffusion of the lipids in the plane of the bilayer are calculated and compared with experimental values. An optimal set of simulation parameters for which the GROMOS96 parameter set 45A3 yields a head group area, chain order parameters and a lateral diffusion coefficient in accordance with the experimental data is listed.

Photoactive yellow protein (PYP) is a prototype of the PAS domain superfamily of signaling proteins. The signaling process is coupled to a three-state photocycle. After the photoinduced trans-cis isomerization of the chromophore, 4-hydroxycinnamic acid (pCA), an early intermediate (pR) is formed, which proceeds to a second intermediate state (pB) on a sub-millisecond time scale. The signaling process is thought to be connected to the conformational changes upon the formation of pB and its recovery to the ground state (pG), but the exact signaling mechanism is not known. Experimental studies of PYP by solution NMR and X-ray crystallography suggest a very flexible protein backbone in the ground as well as in the signaling state. The relaxation from the pR to the pB state is accompanied by the protonation of the chromophore's phenoxyl group. This was found to be of crucial importance for the relaxation process. With the goal of gaining a better understanding of these experimental observations on an atomistic level, we performed five MD simulations on the three different states of PYP: a 1 ns simulation of PYP in its ground state [pG(MD)], a 1 ns simulation of the pR state [pR(MD)], a 2 ns simulation of the pR state with the chromophore protonated (pRprot), a 2 ns simulation of the pR state with Glu46 exchanged by Gln (pRGln) and a 2 ns simulation of PYP in its signaling state [pB(MD)]. Comparison of the pG simulation results with X-ray and NMR data, and with the results obtained for the pB simulation, confirmed the experimental observations of a rather flexible protein backbone and conformational changes during the recovery of the pG from the pB state. The conformational changes in the region around the chromophore pocket in the pR state were found to be crucially dependent on the strength of the Glu46-pCA hydrogen bond, which restricts the mobility of the chromophore in its unprotonated form considerably. Both the mutation of Glu46 with Gln and the protonation of the chromophore weaken this hydrogen bond, leading to an increased mobility of pCA and large structural changes in its surroundings. These changes, however, differ considerably during the pRGln and pRprot simulations, providing an atomistic explanation for the enhancement of the rate constant in the Gln46 mutant.

A comparison of six GROMOS parameter sets for aliphatic alkanes is presented using simulations of a triglyceride, trioctanoin, in the gel or alpha phase. It is found that the parameter set 43A2 results in a system that is rigid and close packed, forcing a collective tilt of the molecules, known experimentally to be disallowed. The parameter set 45A3-45 results in an expanded system where the molecules undergo a conformational change from the chair to tuning-fork form. The lamellae are fragmented, with the fragments shifted in the direction parallel to the layer normal. The newest parameter set 45A3, when used with either of two sets of partial charges for the glycerol ester group, retains the characteristics of the gel phase. The last parameter set tested, 45A3-45x12, also performs well when combined with either set of partial charges, although the system expands in the plane of the bilayer.

Aggregates of n-dodecyl phosphate present an attractive model system of simple phospholipid amphiphile supramolecular structures for study by molecular dynamics simulation, since these systems have previously been studied experimentally under various conditions. A detailed molecular dynamics description of the properties of planar bilayer membranes (as a model for unilamellar vesicular membranes) and spherical micelles under various simulated conditions is presented. It is shown that the united-atom model of GROMOS96 applying the force-field parameter set 43A2 for biomolecular systems yields properties in agreement with experimental ones in most cases. Hydrogen bonding plays a role in stabilizing the bilayer aggregates at low pH, but not for the micelles, which are energetically favoured at high pH. NMR -S(CD) order parameters for a lipid bilayer system, the diffusion of amphiphiles within aggregates and of counterions, and lifetimes of hydrogen bonds between amphiphiles and to water are estimated from the MD simulations.

Insulin crystallizes in different forms, some of which show different conformations for the different molecules in the asymmetric unit. This observation leads to the question as to which conformation the molecule will adopt in solution. Molecular dynamics computer simulations of rhombohedral 2 Zn pig insulin have been carried out for both monomers (1 and 2) independently in order to study their behaviour in the absence of quaternary structure and crystal packing forces. These preliminary 120 ps simulations suggest that both monomers converge in solution to very similar conformations which differ from the X-ray structures of both monomer 1 and 2 (Chinese nomenclature), but are closer to the former, as has previously been suggested by an analysis of the crystal packing (Chothia et al. 1983) and by energy minimization (Wodak et al. 1984). The secondary structure of the molecules is basically preserved, as expected. A detailed description of the conformational changes is given.

The structure of crystalline alpha-cyclodextrin (alpha-CD) hexahydrate, form I (C36H60O30.6H2O, space group P2(1)2(1)2(1)) is experimentally so well determined by X-ray and by neutron diffraction analyses that the positions of all the hydrogen atoms are available. This provides an opportunity for testing an empirical force field that is currently used in simulations of protein and nucleic acid structures by performing molecular dynamics studies employing the GROMOS program package on a system of 4 unit cells containing 16 alpha-CD molecules and 96 water molecules. A detailed comparison of the simulated and experimentally determined crystal structures shows that the experimental positions of the alpha-CD atoms are reproduced within 0.025 nm, well within the overall experimental accuracy of 0.036 nm; that the water molecules are on average within 0.072 nm from their experimental sites, with two thirds reproduced within experimental accuracy by the calculations; that high correlation is produced between the occurrence of simulated and experimentally observed hydrogen bonds. The good agreement between simulated and experimental results suggests that the tested force field is reliable.

Molecular dynamics (MD) simulations for crystalline beta-cyclodextrin dodecahydrate (beta-CD) at two different temperatures, 293 K and 120 K, have been performed using the GROMOS program package. The calculated structural properties are compared to those obtained from neutron diffraction studies of this system at the quoted temperatures. The simulation was carried out over a period of 20 ps on four unit cells containing 8 beta-CD molecules and 96 water molecules, whereby all atoms were allowed to move. At room temperature, the experimental positions of the (non-hydrogen) glucose atoms are reproduced within 0.034 nm, a value which is smaller than the experimental (0.041 nm) or simulated (0.049 nm) overall root mean square (rms) positional fluctuation. The corresponding numbers for the low temperature study are 0.046 nm, 0.019 nm and 0.022 nm. At both temperatures the experimentally observed degree of anisotropy of the atomic motions is also found in the simulations. The comparison of a variety of structural properties leads to the conclusion that the molecular model and force field used are able to simulate the cyclodextrin system very well. Experimentally observed differences in properties as a function of number of glucose units in the CD molecule (alpha-CD, 6 versus beta-CD, 7) and as a function of temperature are qualitatively reproduced by the simulations.

A molecular dynamics simulation was carried out with avian pancreatic polypeptide hormone (aPP) as an isolated monomer explicitly including the solvent (MDS). The simulation and the resulting mean structure are compared with the results of a corresponding crystal simulation (MDC) with 4 aPP molecules plus interstitial water in a periodic boundary unit cell and with the X-ray structure (van Gunsteren, Haneef et al., manuscript in preparation). Comparison is based on the time span 5 to 15 ps and considering cartesian coordinates, dihedral angles, H-bond length, and accessible surface area. While in the MDC simulation equilibration is fast and complete, it does occur in MDS for most but not all parts of the molecule; the turn region starts moving away from the X-ray structure after 9 ps. Only minor differences result when dimer-forming side chains, e.g. tyrosines 7 and 21, are exposed to solvent. The largest rms fluctuations are encountered in exposed polar side chains of Asp 11, Glu 15, Arg 19, and Arg 33, but also in the hydrophobic core residue Phe 20, the only phenylalanine residue present. The latter undergoes an abrupt reorientation suitable for verification by NMR spectroscopy, which is possibly related to the motion of the turn region. The main-chain dihedral angles of the alpha-helix are shifted from values generally found in crystal structures towards those of the ideal Pauling helix. There is concomitant H-bond elongation. The effects are most pronounced and consistent in the MDS simulation.