Search

A method for structure refinement of molecules based on residual dipolar coupling (RDC) data is proposed. It calculates RDC values using rotational and molecule-internal configurational sampling instead of the common refinement procedure that is based on the approximation of the nonuniform rotational distribution of the molecule by a single alignment tensor representing the average nonuniformity of this distribution. Using rotational sampling, as is occurring in the experiment leading to observable RDCs, the method stays close to the experiment. It avoids the use of an alignment tensor and thus the assumption that the overall rotation of the molecule is decoupled from its internal motions and that the molecule be rigid. Two simple molecules, two-united-atomic ethane and a cyclooctane molecule with eight side chains, containing 24 united atoms, serve as the so-called "toy model" test systems. The method demonstrates the influence of molecular flexibility and force-field deficiencies on the outcome of structure refinement based on RDCs. For a molecule of a given size (number of atoms N at ), there must be a sufficiently large number N RDC of measured RDC values available to allow the restraining forces to bias the overall orientation distribution of the molecule. If the ratio N RDC / N at gets too small, the RDC-restraining forces will either not be strong enough to change the overall rotational direction of the molecule such that the target RDC values are approximated well or will be so strong that they induce a local deformation of the molecule. In the latter case, the size or inertia of the molecule hinders a restraining-induced overall rotation and the internal structure of the molecule is not strong enough to avoid local deformation due to the restraining forces.

The derivation of protein structure from values of observable quantities measured in NMR experiments is a rather nontrivial task due to (i) the limited number of data compared to degrees of freedom of a protein, (ii) the uncertainty inherent to the function connecting an observable quantity to molecular structure, (iii) the finite quality of biomolecular models and force fields used in structure refinement, and (iv) the conformational freedom of a protein in aqueous solution, which requires extensive conformational sampling and appropriate conformational averaging when calculating or restraining to sets of NMR data. The protein interleukin-4 (IL-4) has been taken as a test case using NOE distances, S 2 order parameters, and 3 J-couplings as test data and the former two types of data as restraints. It is shown that, by combining sets of different, complementary NMR data as restraints in MD simulations, inconsistencies in the data or flaws in the model and procedures used to derive protein structure from NMR data can be detected. This leads to an improved structural interpretation of such data particularly in more mobile loop regions.

The role of entropic contributions in processes involving biomolecules is illustrated using the process of vaporization or condensation of the solvents water and methanol and the process of polypeptide folding in solution using molecular models at different levels of resolution: subatomic, atomic, supra-atomic, and supramolecular. For the folding process, a β-hexapeptide that adopts, as inferred from NMR experiments, both a right-handed 2.710/12-helical fold and a left-handed 314-helical fold in methanol, is used to illustrate the challenge of modeling thermodynamically driven processes at different levels of resolution.

Starting from a nonpolarizable rigid molecular model for DMSO, a polarizable rigid model for DMSO in the liquid phase is proposed. The molecular polarizability is represented by a charge-on-spring (COS) inducible dipole at the sulfur atom and its polarization is damped for large values of the electric field strength. Some parameters of the model, the partial charge distribution, the polarizability, the value of the electric field at which damping sets in, were varied, and the overall depth and repulsion strength of the Lennard-Jones interactions were scaled with the aim of reproducing the experimental values for the molecular dipole moment, the density, heat of vaporization, and static dielectric constant of liquid DMSO at ambient temperature and pressure. The polarizable DMSO/D model well reproduces the experimental data, such as excess free energy, surface tension, translational, and rotational diffusion, for liquid DMSO. The model is further tested in mixtures with polarizable water models and in mixtures with nonpolarizable DMSO. The latter results suggest that it could serve as model for polarizable cosolvent in aqueous solutions.

The predictive power of biomolecular simulation critically depends on the quality of the force field or molecular model used and on the extent of conformational sampling that can be achieved. Both issues are addressed. First, it is shown that widely used force fields for simulation of proteins in aqueous solution appear to have rather different propensities to stabilize or destabilize α-, π-, and 3(10)- helical structures, which is an important feature of a biomolecular force field due to the omni-presence of such secondary structure in proteins. Second, the relative stability of secondary structure elements in proteins can only be computationally determined through so-called free-energy calculations, the accuracy of which critically depends on the extent of configurational sampling. It is shown that the method of enveloping distribution sampling is a very efficient method to extensively sample different parts of configurational space.

NMR experiments have shown that water molecules in the crystal of the protein Crh are still mobile at temperatures well below 273 K. In order to investigate this water anomaly, a molecular dynamics (MD) simulation study of crystalline Crh was carried out to determine the mobility of water in this crystal. The simulations were carried out at three temperatures, 150, 200, and 291 K. Simulations of bulk water at these temperatures were also done to obtain the properties of the simple point charge (SPC) water model used at these temperatures and to allow a comparison of the properties of water in the Crh crystal with those of bulk water at the same temperatures. According to the simulations, water is immobilized at 150 K both in crystal and in bulk water. As expected, at 291 K it diffuses and rotates more slowly in the protein crystal than in bulk water. However, at 200 K, the translational and rotational mobility of the water molecules is larger in the crystal than in bulk water. The enhancement of water mobility in the crystal at 200 K was further investigated by MD simulations in which the backbone or all protein atoms were positionally restrained, and in which additionally the electrostatic protein-water interactions were removed. Of these changes in the environment of the water molecules, rigidifying the protein backbones slightly enhanced water diffusion, while it slowed down rotation. In contrast, removal of electrostatic protein-water interactions did not change water diffusion but enhanced rotational motion significantly. Further investigations are required to delineate particular features of the protein crystal that induce the anomalous behavior of water at 200 K.

Atomistic molecular dynamics simulations of proteins in aqueous solution are still limited to the multinanosecond time scale and multinanometer range by computational cost. Combining atomic solutes with a supra-molecular solvent model in hybrid fine-grained/coarse-grained (FG/CG) simulations allows atomic detail in the region of interest while being computationally more efficient. A recent comparison of the properties of four proteins in CG water versus FG water showed the preservation of the secondary and tertiary structure with a computational speed-up of at least an order of magnitude. However, an increased occurrence of hydrogen bonds between side chains was observed due to a lack of hydrogen-bonding partners in the supra-molecular solvent. Here, the introduction of a FG water layer around the protein to recover the hydrogen-bonding pattern of the atomistic simulations is studied. Three layer thicknesses of 0.2, 0.4, and 0.8 nm are considered. A layer thickness of 0.8 nm is found sufficient to recover the behavior of the proteins in the atomistic simulations, whereas the hybrid simulation is still three times more efficient than the atomistic one and the cutoff radius for nonbonded interactions could be increased from 1.4 to 2.0 nm.

The calculation of the relative free energies of ligand-protein binding, of solvation for different compounds, and of different conformational states of a polypeptide is of considerable interest in the design or selection of potential enzyme inhibitors. Since such processes in aqueous solution generally comprise energetic and entropic contributions from many molecular configurations, adequate sampling of the relevant parts of configurational space is required and can be achieved through molecular dynamics simulations. Various techniques to obtain converged ensemble averages and their implementation in the GROMOS software for biomolecular simulation are discussed, and examples of their application to biomolecules in aqueous solution are given., (© 2011 American Chemical Society)

The ability to design well-folding β-peptides with a specific biological activity requires detailed insight into the relationship between the β-amino acid sequence and the dominant three-dimensional structure of such a peptide. To this end, secondary structure preferences of two sets of 16 β-peptides were investigated by means of one-step perturbation using molecular dynamics (MD) simulations. For each set of peptides, two reference-state simulations and one perturbed-state simulation were carried out to predict the secondary structure preferences for the other 15 peptides. The results show that the substitution of a methyl group in the third or fourth residue stabilizes the left-handed 3(14)-helix over the right-handed 2.7(10/12)-helix for the set of hexapeptides A; for the set of heptapeptides B, having methyl substitutions at both β- and α-carbon positions of the fourth or fifth residue stabilizes the left-handed 3(14)-helix over the right-handed 2.5(12)-helix. Not only the side-chain substitution pattern but also the side-chain composition affects the relative stability of different secondary structures. The approach described here may be of use in peptide design with an eye to obtaining peptides with particular folds and biological activities.

A refined implicit aqueous solvation model is proposed for the simulation of biomolecules without the explicit inclusion of the solvent degrees of freedom. The mean force due to solvation is approximated by the derivative of a simple analytic function of the solvent accessible surface area combined with two atomic solvation parameters, as described previously, with the addition of a novel term to account for the interaction of the interior atoms of the solute with the solvent. The extended model is parametrized by comparing the structural properties and energies computed from simulations of six test proteins of varying sizes and shapes using the new solvation energy term with the corresponding values obtained from simulations in vacuum, using the original implicit solvent model and in explicit water, and from the X-ray or NMR model structures. The mean solvation model proposed here improves the structural properties relative to vacuum simulations and relative to the simpler model that neglects the volume contribution, while remaining significantly more efficient than simulations in explicit water.

Three methods to enhance the configurational sampling of ions in aqueous solution, temperature and Hamiltonian replica exchange and adiabatic decoupling with translational temperature scaling, were compared for a system of CaSO(4) in water. It took 11 replicas in the case of temperature replica exchange to make use of a diffusion coefficient that is a factor of 1.5 larger at 350 K compared to that at 300 K. Thirty replicas were required in the Hamiltonian replica exchange with charge reduction to reach uncharged ions that have a diffusion coefficient that is 2-7 times larger than the fully charged ones. The adiabatic decoupling technique with translational temperature scaling yielded a diffusion coefficient that was 15 times larger while keeping the distribution of the water molecules around the ions unaltered with respect to the standard temperature simulation. This result illustrates the efficiency of the adiabatic decoupling technique to enhance configurational sampling.

With the aim to gain a better understanding of the various driving forces that govern sequence specific DNA minor groove binding, we performed a thermodynamic analysis of netropsin binding to an AT-containing and to a set of six mixed AT/GC-containing binding sequences in the DNA minor groove. The relative binding free energies obtained using molecular dynamics simulations and free energy calculations show significant variations with the binding sequence. While the introduction of a GC base pair in the middle or close to the middle of the binding site is unfavorable for netropsin binding, a GC base pair at the end of the binding site appears to have no negative influence on the binding. The results of the structural and energetic analyses of the netropsin-DNA complexes reveal that the differences in the calculated binding affinities cannot be explained solely in terms of netropsin-DNA hydrogen-bonding or interaction energies. In addition, solvation effects and entropic contributions to the relative binding free energy provide a more complete picture of the various factors determining binding. Analysis of the relative binding entropy indicates that its magnitude is highly sequence-dependent, with the ratio |TDeltaDeltaS|/|DeltaDeltaH| ranging from 0.07 for the AAAGA to 1.7 for the AAGAG binding sequence, respectively.

Explicit solvent molecular dynamics (MD) simulations of three beta-heptapeptides with a central beta- HAla(alpha-F) amino acid (Figure 1) in methanol are reported. They aim at an analysis of the conformational consequences of C(alpha) carbon atom bound fluoro atoms, and the particular configuration of the central fluoro-beta-amino acid: peptide 3 with an S configuration of the C(alpha) bound fluor atom, peptide 4 with an R configuration of the C(alpha) bound fluor atom, and peptide 5 with a difluoro substitution at the C(alpha) atom of residue 4. The NMR and CD spectra of these three beta-peptides were earlier (Mathad et al. Helv. Chim. Acta 2005, 88, 266-280) interpreted to indicate a decrease in propensity of 3(14)-helical structure from peptide 4 to peptide 5 to peptide 3. This result was at odds with previous experimental data for beta-heptapeptides with a central beta-HAla(alpha-Me) amino acid which showed that the beta-heptapeptide with the S,S configuration of the central beta-HAla(alpha-Me) was the most 3(14)-helical, whereas the S,R configuration did not lead to any detected helicity. The reported MD simulations resolve this paradox. The MD trajectories of all three peptides do agree with the primary, measured data: NMR nuclear Overhauser effect (NOE) atom-atom distance bounds and (3)J-coupling constants. A conformational analysis of the MD trajectory conformations shows, however, a decrease in 3(14)-helical character from peptide 3 to peptide 5 to peptide 4, which is in line with the results for the nonfluorinated peptides. It is shown that interpretation of NMR NOE data using single-structure refinement in vacuo based on local (along the sequence) and limited atom-atom distance data as in ref 1 (Mathad et al. Helv. Chim. Acta 2005, 88, 266-280) may lead to molecular structures that are not representative for the ensemble of molecular conformations.

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.

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.

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.

This article investigates the convergence of structural and dynamical properties with system size and with time in molecular dynamics simulations of solvated phospholipid bilayers performed at constant volume under periodic boundary conditions using lattice-sum electrostatics. The electron density profile across the bilayer, the carbon-deuterium order parameters, and the surface tension are shown to be converged for a bilayer containing 36 lipids per leaflet and simulated over a period of 3-4 ns. Reasonable estimates for these properties can already be obtained from a system containing 16 lipids per leaflet. The convergence limit of 36 lipids per leaflet and the investigation of the correlation between lipid headgroup dipoles suggest a correlation length of about 3-5 nm in the lateral directions for a hydrated DPPC bilayer in the liquid-crystalline phase. Although these (relatively small) system sizes and (relatively short) time scales appear sufficient to obtain converged collective structural properties at constant volume, two restrictions should be kept in mind: (i) the relaxation times associated with the motion of individual lipids may be much longer and (ii) simulated properties converge significantly faster under constant volume conditions as compared to constant pressure conditions. Therefore, an accurate assessment of the dynamical properties of the system or of the relaxation of the bilayer under constant pressure conditions may require longer simulation time scales.