Search

More than a half century ago it became feasible to simulate, using classical-mechanical equations of motion, the dynamics of molecular systems on a computer. Since then classical-physical molecular simulation has become an integral part of chemical research. It is widely applied in a variety of branches of chemistry and has significantly contributed to the development of chemical knowledge. It offers understanding and interpretation of experimental results, semiquantitative predictions for measurable and nonmeasurable properties of substances, and allows the calculation of properties of molecular systems under conditions that are experimentally inaccessible. Yet, molecular simulation is built on a number of assumptions, approximations, and simplifications which limit its range of applicability and its accuracy. These concern the potential-energy function used, adequate sampling of the vast statistical-mechanical configurational space of a molecular system and the methods used to compute particular properties of chemical systems from statistical-mechanical ensembles. During the past half century various methodological ideas to improve the efficiency and accuracy of classical-physical molecular simulation have been proposed, investigated, evaluated, implemented in general simulation software or were abandoned. The latter because of fundamental flaws or, while being physically sound, computational inefficiency. Some of these methodological ideas are briefly reviewed and the most effective methods are highlighted. Limitations of classical-physical simulation are discussed and perspectives are sketched.

A method for structure refinement of molecules based on residual dipolar coupling (RDC) data is proposed. It calculates RDC values using magnetic-field rotational sampling of the rotational degrees of freedom of a molecule in conjunction with molecule-internal configurational sampling. By applying 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, therefore, the assumptions that the overall rotation of the molecule is decoupled from its internal motions and that the molecule is rigid. Two simple molecules, a relatively rigid and a very flexible cyclo-octane molecule with eight aliphatic side chains containing 24 united atoms, serve as so-called "toy model" test systems. The method demonstrates the influence of molecular flexibility, force-field dominance, and the number of RDC restraints available on the outcome of structure refinement based on RDCs. Magnetic-field rotational sampling is basically equivalent but more efficient than explicitly sampling the rotational degrees of freedom of the molecule. In addition, the performance of the method is less dependent on the number NRDC of measured RDC-values available. The restraining forces bias the overall orientation distribution of the molecule correctly. This study suggests that the information content of RDCs with respect to molecular structure is limited., (© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/).)

In NMR experiments, residual dipolar couplings (RDCs) in a molecule can be measured by averaging the dipolar couplings (DCs) over the rotational motion of a molecule in an environment that induces a slight anisotropic orientation distribution of the molecule. Since the shape of the anisotropic distribution cannot be measured, it is standard practice to use a particular orientation distribution of the molecule with respect to the magnetic field, in the form of a so-called alignment tensor (AT), to calculate RDC-values for the molecule. Since the same alignment tensor is commonly used to calculate the different RDCs of a molecule, this approach rests on the assumption that the rotational motion of the molecule is decoupled from its internal motions and that the molecule is rigid. The validity of these two assumptions is investigated for a small, simple molecule, using a relatively rigid atomic interaction function or force field and a more flexible one. By simulating the molecule using an orientation-biasing force an anisotropic rotational distribution can be generated, for which RDCs can be obtained. Using these RDCs as target RDCs when applying one of the two approaches of structure refinement based on RDCs, it can be investigated how well the target RDCs are approximated in the RDC restraining and whether the corresponding nonuniform orientation distribution is reproduced. For the relatively rigid version of the molecule, the AT approach reproduces the target RDC-values, although the nonuniform orientation distribution of the angle θ ab,H between the vector r⃗ ab connecting two atoms a and b in the molecule and the vector representing the direction of the magnetic field H⃗ as generated in the orientation-biasing simulation cannot be reproduced in the AT RDC-restraining simulation. For the relatively flexible version of the molecule, the AT approach fails to reproduce both the target RDC values and the nonuniform orientation distribution. For biomolecules with flexible parts, the application of the AT approach is thus not recommended. Instead, a method based on sampling of the rotational and internal degrees of freedom of the molecule should be applied in molecular structure determination or refinement based on measured RDCs.

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.

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).)

Computer simulation of proteins in aqueous solution at the atomic level of resolution is still limited in time span and system size due to limited computing power available and thus employs a variety of time-saving techniques that trade some accuracy against computational effort. An example of such a time-saving technique is the application of constraints to particular degrees of freedom when integrating Newton's or Langevin's equations of motion in molecular dynamics (MD) or stochastic dynamics (SD) simulations, respectively. The application of bond-length constraints is standard practice in protein simulations and allows for a lengthening of the time step by a factor of three. Applying recently proposed algorithms to constrain bond angles or dihedral angles, it is investigated, using the protein trypsin inhibitor as test molecule, whether bond angles and dihedral angles involving hydrogen atoms or even stiff proper (torsional) dihedral angles as well as improper ones (maintaining particular tetrahedral or planar geometries) may be constrained without generating too many artificial side effects. Constraining the relative positions of the hydrogen atoms in the protein allows for a lengthening of the time step by a factor of two. Additionally constraining the improper dihedral angles and the stiff proper (torsional) dihedral angles in the protein does not allow for an increase of the MD or SD time step., (© 2021 The Authors. Proteins: Structure, Function, and Bioinformatics published by Wiley Periodicals LLC.)

Computer simulation of proteins in aqueous solution at the atomic level of resolution is still limited in time span and system size due to limited computing power available and thus employs a variety of time-saving techniques that trade some accuracy against computational effort. Examples of such time-saving techniques are the application of constraints to particular degrees of freedom or the use of a multiple-time-step (MTS) algorithm distinguishing between particular forces when integrating Newton's equations of motion. The application of two types of MTS algorithms to bond-stretching forces versus the remaining forces in molecular dynamics (MD) simulations of a protein in aqueous solution or of liquid water is investigated and the results in terms of total energy conservation and the influence on various other properties are compared to those of MD simulations of the same systems using bond-length, and for water bond-angle, constraints. At comparable computational effort, the use of bond-length constraints in proteins leads to better energy conservation and less distorted properties than the two MTS algorithms investigated., (© 2021 The Authors. Journal of Computational Chemistry published by Wiley Periodicals LLC.)

Values of S 2 CH and S 2 NH order parameters derived from NMR relaxation measurements on proteins cannot be used straightforwardly to determine protein structure because they cannot be related to a single protein structure, but are defined in terms of an average over a conformational ensemble. Molecular dynamics simulation can generate a conformational ensemble and thus can be used to restrain S 2 CH and S 2 NH order parameters towards experimentally derived target values S 2 CH (exp) and S 2 NH (exp). Application of S 2 CH and S 2 NH order-parameter restraining MD simulation to bond vectors in 63 side chains of the protein hen egg white lysozyme using 51 S 2 CH (exp) target values and 28 S 2 NH (exp) target values shows that a conformational ensemble compatible with the experimentally derived data can be obtained by using this technique. It is observed that S 2 CH order-parameter restraining of C-H bonds in methyl groups is less reliable than S 2 NH order-parameter restraining because of the possibly less valid assumptions and approximations used to derive experimental S 2 CH (exp) values from NMR relaxation measurements and the necessity to adopt the assumption of uniform rotational motion of methyl C-H bonds around their symmetry axis and of the independence of these motions from each other. The restrained simulations demonstrate that side chains on the protein surface are highly dynamic. Any hydrogen bonds they form and that appear in any of four different crystal structures, are fluctuating with short lifetimes in solution., (© 2020 The Authors. ChemBioChem published by Wiley-VCH GmbH.)

An algorithm to apply bond-angle constraints in molecular dynamics simulations of macromolecules or molecular liquids is presented. It uses Cartesian coordinates and determines the Lagrange multipliers required for maintaining the constraints iteratively. It constitutes an alternative to the use of only distance constraints (DCs) between particles to maintain a particular geometry. DCs are unsuitable to maintain particular, for example, linear or flat, geometries of molecules. The proposed algorithm can easily handle bond-length, bond-angle, and dihedral-angle constraints simultaneously, as when calculating a potential of mean force along a dihedral-angle degree of freedom., (© 2020 Wiley Periodicals LLC.)

Computer simulations of molecular systems enable structure-energy-function relationships of molecular processes to be described at the sub-atomic, atomic, supra-atomic or supra-molecular level and plays an increasingly important role in chemistry, biology and physics. To interpret the results of such simulations appropriately, the degree of uncertainty and potential errors affecting the calculated properties must be considered. Uncertainty and errors arise from (1) assumptions underlying the molecular model, force field and simulation algorithms, (2) approximations implicit in the interatomic interaction function (force field), or when integrating the equations of motion, (3) the chosen values of the parameters that determine the accuracy of the approximations used, and (4) the nature of the system and the property of interest. In this overview, advantages and shortcomings of assumptions and approximations commonly used when simulating bio-molecular systems are considered. What the developers of bio-molecular force fields and simulation software can do to facilitate and broaden research involving bio-molecular simulations is also discussed., (© 2020 Wiley-VCH GmbH.)

Values of 3 J-couplings as obtained from NMR experiments on proteins cannot easily be used to determine protein structure due to the difficulty of accounting for the high sensitivity of intermediate 3 J-coupling values (4-8 Hz) to the averaging period that must cover the conformational variability of the torsional angle related to the 3 J-coupling, and due to the difficulty of handling the multiple-valued character of the inverse Karplus relation between torsional angle and 3 J-coupling. Both problems can be solved by using 3 J-coupling time-averaging local-elevation restraining MD simulation. Application to the protein hen egg white lysozyme using 213 backbone and side-chain 3 J-coupling restraints shows that a conformational ensemble compatible with the experimental data can be obtained using this technique, and that accounting for averaging and the ability of the algorithm to escape from local minima for the torsional angle induced by the Karplus relation, are essential for a comprehensive use of 3 J-coupling data in protein structure determination.

Various algorithms to apply dihedral-angle constraints in molecular dynamics or stochastic dynamics simulations of molecular systems are presented, investigated, and tested. They use Cartesian coordinates and determine the Lagrangian multipliers necessary for maintaining the constraints iteratively. The most suitable algorithm to maintain a dihedral-angle constraint is numerically compared to the alternative to use distance constraints to this end. It can easily be used to obtain a potential of mean force along a dihedral-angle coordinate.

Epothilones are among the most potent chemotherapeutic drugs used for the treatment of cancer. Epothilone A (EpoA), a natural product, is a macrocyclic molecule containing 34 non-hydrogen atoms and a thiazole side chain. NMR studies of EpoA in aqueous solution, unbound as well as bound to αβ-tubulin, and unbound in dimethyl sulfoxide (DMSO) solution have delivered sets of nuclear Overhauser effect (NOE) atom-atom distance bounds, but no structures based on NMR data are present in structural data banks. X-ray diffraction of crystals has provided structures of EpoA unbound and bound to αβ-tubulin. Since both crystal structures derived from X-ray diffraction intensities do not completely satisfy the three available sets of NOE distance bounds for EpoA, molecular dynamics (MD) simulations have been employed to obtain conformational ensembles in aqueous and in DMSO solution that are compatible with the respective NOE data. It was found that EpoA displays a larger conformational variability in DMSO than in water and the two conformational ensembles show little overlap. Yet, they both provide conformational scaffolds that are energetically accessible at physiological temperature and pressure.

Computer simulation of molecular systems enables structure-energy-function relationships of molecular processes to be described at the sub-atomic, atomic, supra-atomic, or supra-molecular level. To interpret results of such simulations appropriately, the quality of the calculated properties must be evaluated. This depends on the way the simulations are performed and on the way they are validated by comparison to values Q exp of experimentally observable quantities Q. One must consider 1) the accuracy of Q exp , 2) the accuracy of the function Q(r N ) used to calculate a Q-value based on a molecular configuration r N of N particles, 3) the sensitivity of the function Q(r N ) to the configuration r N , 4) the relative time scales of the simulation and experiment, 5) the degree to which the calculated and experimental properties are equivalent, and 6) the degree to which the system simulated matches the experimental conditions. Experimental data is limited in scope and generally corresponds to averages over both time and space. A critical analysis of the various factors influencing the apparent degree of (dis)agreement between simulations and experiment is presented and illustrated using examples from the literature. What can be done to enhance the validation of molecular simulation is also discussed., (© 2018 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim.)

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.

At low pH, human growth hormone (hGH) adopts a partially folded state, in which the native helices are maintained, but the long loop regions and side-chain packing become disordered. Some of the S 2 order parameters for backbone N-H vectors derived from NMR relaxation measurements on hGH at low pH initially seem contradictory. Three isolated residues (15, 20, and 171) in helices A and D exhibit low order parameter values (<0.5) indicating flexibility, whereas residue 143 in the centre of a long flexible loop region has a high order parameter (0.82). Using S 2 order parameter restraining MD simulations, this paradox has been resolved. Low S 2 values in helices are due to the presence of a mixture of 3 10 -helical and α-helical hydrogen bonds. High S 2 values in relatively disordered parts of a protein may be due to fluctuating networks of hydrogen bonds between the backbone and the side chains, which restrict the motion of N-H bond vectors., (© 2017 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim.)

During the past half century, the number and accuracy of experimental techniques that can deliver values of observables for biomolecular systems have been steadily increasing. The conversion of a measured value Q exp of an observable quantity Q into structural information is, however, a task beset with theoretical and practical problems: 1) insufficient or inaccurate values of Q exp , 2) inaccuracies in the function Q(r→) used to relate the quantity Q to structure r→ , 3) how to account for the averaging inherent in the measurement of Q exp , 4) how to handle the possible multiple-valuedness of the inverse r→(Q) of the function Q(r→) , to mention a few. These apply to a variety of observable quantities Q and measurement techniques such as X-ray and neutron diffraction, small-angle and wide-angle X-ray scattering, free-electron laser imaging, cryo-electron microscopy, nuclear magnetic resonance, electron paramagnetic resonance, infrared and Raman spectroscopy, circular dichroism, Förster resonance energy transfer, atomic force microscopy and ion-mobility mass spectrometry. The process of deriving structural information from measured data is reviewed with an eye to non-experts and newcomers in the field using examples from the literature of the effect of the various choices and approximations involved in the process. A list of choices to be avoided is provided., (© 2016 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim.)

The cyclic octa-peptide octreotide and its derivatives are used as diagnostics and therapeutics in relation to particular types of cancers. This led to investigations of their conformational properties using spectroscopic, NMR and CD, methods. A CF 3 -substituted derivative, that was designed to stabilize the dominant octreotide conformer responsible for receptor binding, turned out to have a lower affinity. The obtained spectroscopic data were interpreted as to show an increased flexibility of the CF 3 derivative compared to the unsubstituted octreotide, which could then explain the lower affinity. In this article, we use MD simulation without and with time-averaged NOE distance and time-averaged local-elevation 3 J-coupling restraining representing experimental NMR data to determine the conformational properties of the different peptides in the different solvents for which experimental data are available, that are compatible with the NOE atom-atom distance bounds and the 3 J HNHα -couplings as derived from the NMR measurements. The conformational ensembles show that the CF 3 substitution in combination with the change of solvent from water to methanol leads to a decrease in flexibility and a shift in the populations of the dominant conformers that are compatible with the experimental data., (Copyright © 2016 Elsevier Ltd. All rights reserved.)

Deriving molecular structure from [Formula: see text]-couplings obtained from NMR experiments is a challenge due to (1) the uncertainty in the Karplus relation [Formula: see text] connecting a [Formula: see text]-coupling value to a torsional angle [Formula: see text], (2) the need to account for the averaging inherent to the measurement of [Formula: see text]-couplings, and (3) the sampling road blocks that may emerge due to the multiple-valuedness of the inverse function [Formula: see text] of the function [Formula: see text]. Ways to properly handle these issues in structure refinement of biomolecules are discussed and illustrated using the protein hen egg white lysozyme as example.

The structural variability of a 16-residue loop (residues 246-261) which is in part disordered and connects two layers of the β-solenoid formed by the prion-form of HET-s and its prion domain HET-s(218-289) is investigated using molecular dynamics computer simulation. A system of three HET-s(218-289) molecules in a β-sheet structure as in the fibril is simulated in aqueous solution. The trajectory structures appear to be consistent with the Cα chemical shift data obtained. In order to delineate the influence of the β-sheet core of the fibril upon the structural variability of the loop, the latter is also simulated without the β-sheet core, but with its N- and C-terminal residues restrained at their positions in the fibril. The analysis of the trajectories shows that the structural variability of the loop is restricted by the β-sheet core, least at its N-terminal end and most in the middle of the trimer.

The multistep umbrella sampling method, which belongs to pathway-dependent methods to calculate conformational free enthalpy differences, is used to calculate the free enthalpy difference between a right-handed 2.710/12 -helix and a left-handed 314 -helix of a hexa-β-peptide in methanol solution. The same conformational free enthalpy difference was previously investigated using pathway-independent methods such as direct counting and enveloping distribution sampling. Our results show that the pathway-dependent simulations are sensitive to the choice of the pathway and its parameter values. A pathway based on restraining distances of hydrogen-bonding atom pairs shows poor sampling for two different values of the restraining force constant. Another pathway based on restraining backbone dihedral angles did smoothly sample the transition between the two helical conformations, but only with a proper choice of the restraining force constant. The results illustrate that if, and only if, a proper pathway and proper parameters are chosen, the multistep umbrella sampling can be almost 50 times more efficient than the pathway-independent methods in this case. The analysis illustrates the advantages and pitfalls of the much used multistep umbrella sampling methodology., (© 2015 The Protein Society.)

An algorithm is proposed for performing molecular dynamics (MD) simulations of a biomolecular solute represented at atomistic resolution surrounded by a surface layer of atomistic fine-grained (FG) solvent molecules within a bulk represented by coarse-grained (CG) solvent beads. The method, called flexible boundaries for multiresolution solvation (FBMS), is based on: (i) a three-region layering of the solvent around the solute, involving an FG layer surrounded by a mixed FG-CG buffer layer, itself surrounded by a bulk CG region; (ii) a definition of the layer boundary that relies on an effective distance to the solute surface and is thus adapted to the shape of the solute as well as adjusts to its conformational changes. The effective surface distance is defined by inverse-nth power averaging over the distances to all non-hydrogen solute atoms, and the layering is enforced by means of half-harmonic distance restraints, attractive for the FG molecules and repulsive for the CG beads. A restraint-free region at intermediate distances enables the formation of the buffer layer, where the FG and CG solvents can mix freely. The algorithm is tested and validated using the GROMOS force field and the associated FG (SPC) and CG (polarizable CGW) water models. The test systems include pure-water systems, where one FG molecule plays the role of a solute, and a deca-alanine peptide with two widely different solute shapes considered, α-helical and fully extended. In particular, as the peptide unfolds, the number of FG molecules required to fill its close-range solvation layer increases, with the additional molecules being provided by the buffer layer. Further validation involves simulations of four proteins in multiresolution FG/CG mixtures. The resulting structural, energetic, and solvation properties are found to be similar to those observed in corresponding pure FG simulations.

In a molecular dynamics (MD) simulation, various thermostat algorithms, including Langevin dynamics (LD), Nosé-Hoover (NH), and weak-coupling (WC) thermostats, can be used to keep the simulation temperature constant. A canonical ensemble is generated by the use of LD and NH, while the nature of the ensemble produced by WC has not yet been identified. A few years ago, it was shown that when using a WC thermostat with particular values of the temperature coupling time for liquid water at ambient temperature and pressure, the distribution of the potential energy is less wide than the canonical one. This led to an artifact in temperature replica-exchange molecular dynamics (T-REMD) simulations in which the potential energy distributions appear not to be equal to the ones of standard MD simulations. In this paper, we re-investigate this problem. We show that this artifact is probably due to the ensemble generated by WC being incompatible with the T-REMD replica-exchange criterion, which assumes a canonical configurational ensemble. We also show, however, that this artifact can be reduced or even eliminated by particular choices of the temperature coupling time of WC and the replica-exchange time period of T-REMD, i.e., when the temperature coupling time is chosen very close to the MD time step or when the exchange time period is chosen large enough. An attempt to develop a T-REMD replica-exchange criterion which is likely to be more compatible with the WC configurational ensemble is reported. Furthermore, an exchange criterion which is compatible with a microcanonical ensemble is used in total energy REMD simulations.

A supra-atomic coarse-grained (CG) force field for liquid n-alkanes is presented. The model was calibrated using experimental thermodynamic data and structural as well as energetic properties for 14 n-alkanes as obtained from atomistic fine-grained (FG) simulations of the corresponding hydrocarbons using the GROMOS 45A3 biomolecular force field. A variation of the nonbonded force-field parameters obtained from mapping the FG interactions onto the CG degrees of freedom to fit the density and heat of vaporization to experimental values turned out to be mandatory for a correct reproduction of these data by the CG model, while the bonded force-field parameters for the CG model could be obtained from a Boltzmann-weighted fit with some variations with respect to the corresponding properties from the FG simulations mapped onto the CG degrees of freedom. The model presents 6 different CG bead types, for bead sizes from 2 to 4 distinguishing between terminal and nonterminal beads within an alkane chain (end or middle). It contains different nonbonded Lennard-Jones parameters for the interaction of CG alkanes with CG water. The CG alkane model was further tested by comparing predictions of the excess free energy, the self-diffusion constant, surface tension, isothermal compressibility, heat capacity, thermal expansion coefficient, and shear viscosity for n-alkanes to experimental values. The CG model offers a thermodynamically calibrated basis for the development of CG models of lipids.

Force field parameters for polarizable coarse-grained (CG) supra-atomic models of liquid cyclohexane are proposed. Two different bead sizes were investigated, one representing two fine-grained (FG) CH(2)r united atoms of the cyclohexane ring, and one representing three FG CH(2)r united atoms. Electronic polarizability is represented by a massless charge-on-spring particle connected to each CG bead. The model parameters were calibrated against the experimental density and heat of vaporization of liquid cyclohexane, and the free energy of cyclohexane hydration. Both models show good agreement with thermodynamic properties of cyclohexane, yet overestimate the self-diffusion. The dielectric properties of the polarizable models agree very well with experiment., (© 2015 Wiley Periodicals, Inc.)

The free enthalpy of changing a nonpolarizable solvent into a polarizable solvent is calculated for an alanine decapeptide solvated in water, methanol, chloroform, or carbon tetrachloride. Introducing polarizability into a water solvent does not change the relative stability between the α-helix and the β-hairpin, while for methanol and chloroform the α-helix is stabilized by about 1 kJ mol(-1) per residue and for carbon tetrachloride by about 2 kJ mol(-1) per residue. These results suggest that the less polar the solvent is, the more the α-helical structure is stabilized with respect to the β-hairpin structure by the use of a polarizable solvent model instead of a nonpolarizable one. This highlights that inclusion of polarizability in models for less polar and nonpolar solvents or protein environments is as important as, if not more important than, including polarizability in models for liquid water.

Background: The contribution of particular hydrogen bonds to the stability of a protein fold can be investigated experimentally as well as computationally by the construction of protein mutants which lack particular hydrogen-bond donors or acceptors with a subsequent determination of their structural stability. However, the comparison of experimental data with computational results is not straightforward. One of the difficulties is related to the representation of the unfolded state conformation., Methods: A series of molecular dynamics simulations of the 34-residue WW domain of protein Pin1 and 20 amide-to-ester mutants started from the X-ray crystal structure and the NMR solution structure are analysed in terms of backbone-backbone hydrogen bonding and differences in free enthalpies of folding in order to provide a structural interpretation of the experimental data available., Results: The contribution of the different β-sheet hydrogen bonds to the relative stability of the mutants with respect to wild type cannot be directly inferred from experimental thermal denaturation temperatures or free enthalpies of chaotrope denaturation for the different mutants, because some β-sheet hydrogen bonds show sizeable variation in occurrence between the different mutants., Conclusions: A proper representation of unfolded state conformations appears to be essential for an adequate description of relative stabilities of protein mutants., General Significance: The simulations may be used to link the structural Boltzmann ensembles to relative free enthalpies of folding between mutants and wild-type protein and show that unfolded conformations have to be treated with a sufficient level of detail in free energy calculations of protein stability. This article is part of a Special Issue entitled Recent developments of molecular dynamics., (Copyright © 2014 Elsevier B.V. All rights reserved.)

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.

Two one-site polarizable urea models, COS/U and COS/U2, based on the charge-on-spring model are proposed. The models are parametrized against thermodynamic properties of urea-water mixtures in combination with the polarizable COS/G2 and COS/D2 models for liquid water, respectively, and have the same functional form of the inter-atomic interaction function and are based on the same parameter calibration procedure and type of experimental data as used to develop the GROMOS biomolecular force field. Thermodynamic, dielectric, and dynamic properties of urea-water mixtures simulated using the polarizable models are closer to experimental data than using the non-polarizable models. The COS/U and COS/U2 models may be used in biomolecular simulations of protein denaturation.

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.

The accuracy of biomolecular simulations depends to some degree on the accuracy of the water model used to solvate the biomolecules. Because many biomolecules such as proteins are electrostatically rather inhomogeneous, containing apolar, polar, and charged moieties or side chains, a water model should be able to represent the polarisation response to a local electrostatic field, while being compatible with the force field used to model the biomolecules or protein. The two polarisable water models, COS/G2 and COS/D, that are compatible with the GROMOS biomolecular force fields leave room for improvement. The COS/G2 model has a slightly too large dielectric permittivity and the COS/D model displays a much too slow dynamics. The proposed COS/D2 model has four interaction sites: only one Lennard-Jones interaction site, the oxygen atom, and three permanent charge sites, the two hydrogens, and one massless off-atom site that also serves as charge-on-spring (COS) polarisable site with a damped or sub-linear dependence of the induced dipole on the electric field strength for large values of the latter. These properties make it a cheap and yet realistic water model for biomolecular solvation.

The flavoenzyme pyranose dehydrogenase (PDH) from the litter decomposing fungus Agaricus meleagris oxidizes many different carbohydrates occurring during lignin degradation. This promiscuous substrate specificity makes PDH a promising catalyst for bioelectrochemical applications. A generalized approach to simulate all 32 possible aldohexopyranoses in the course of one or a few molecular dynamics (MD) simulations is reported. Free energy calculations according to the one-step perturbation (OSP) method revealed the solvation free energies (ΔGsolv) of all 32 aldohexopyranoses in water, which have not yet been reported in the literature. The free energy difference between β- and α-anomers (ΔGβ-α) of all d-stereoisomers in water were compared to experimental values with a good agreement. Moreover, the free-energy differences (ΔG) of the 32 stereoisomers bound to PDH in two different poses were calculated from MD simulations. The relative binding free energies (ΔΔGbind) were calculated and, where available, compared to experimental values, approximated from Km values. The agreement was very good for one of the poses, in which the sugars are positioned in the active site for oxidation at C1 or C2. Distance analysis between hydrogens of the monosaccharide and the reactive N5-atom of the flavin adenine dinucleotide (FAD) revealed that oxidation is possible at HC1 or HC2 for pose A, and at HC3 or HC4 for pose B. Experimentally detected oxidation products could be rationalized for the majority of monosaccharides by combining ΔΔGbind and a reweighted distance analysis. Furthermore, several oxidation products were predicted for sugars that have not yet been tested experimentally, directing further analyses. This study rationalizes the relationship between binding free energies and substrate promiscuity in PDH, providing novel insights for its applicability in bioelectrochemistry. The results suggest that a similar approach could be applied to study promiscuity of other enzymes.

A method is described that allows experimental S(2) order parameters to be enforced as a time-averaged quantity in molecular dynamics simulations. The two parameters that characterize time-averaged restraining, the memory relaxation time and the weight of the restraining potential energy term in the potential energy function used in the simulation, are systematically investigated based on two model systems, a vector with one end restrained in space and a pentapeptide. For the latter it is shown that the backbone N-H order parameter of individual residues can be enforced such that the spatial fluctuations of quantities depending on atomic coordinates are not significantly perturbed. The applicability to realistic systems is illustrated for the B3 domain of protein G in aqueous solution.

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.

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.

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.

Molecular dynamics simulation of biomolecules in solvent using an atomic model for both the biomolecules and the solvent molecules is still computationally rather demanding considering the time scale of the biomolecular motions. The use of a supramolecular coarse-grained (CG) model can speed up the simulation considerably, but it also reduces the accuracy inevitably. Combining an atomic fine-grained (FG) level of modeling for the biomolecules and a supramolecular CG level for the solvent into a hybrid system, the increased computational efficiency may outweigh the loss of accuracy with respect to the biomolecular properties in the hybrid FG/CG simulation. Here, a previously published CG methanol model is reparametrized, and then a 1:1 mixture of FG and CG methanol is used to calibrate the FG-CG interactions using thermodynamic and dielectric screening data for liquid methanol. The FG-CG interaction parameter set is applied in hybrid FG/CG solute/solvent simulations of the folding equilibria of three β-peptides that adopt different folds. The properties of the peptides are compared with those obtained in FG solvent simulations and with experimental NMR data. The comparison shows that the folding equilibria in the pure CG solvent simulations are different from those in the FG solvent simulations because of the lack of hydrogen-bonding partners in the supramolecular CG solvent. Next, we introduced an FG methanol layer around the peptides in CG solvent to recover the hydrogen-bonding pattern of the FG solvent simulations. The result shows that with the FG methanol layer, the folding equilibria of the three β-peptides are very similar to those in the FG solvent simulations, while the computational efficiency is at least 3 times higher and the cutoff radius for nonbonded interactions could be increased from 1.4 to 2.0 nm.

Electronic polarizability is usually treated implicitly in molecular simulations, which may lead to imprecise or even erroneous molecular behavior in spatially electronically inhomogeneous regions of systems such as proteins, membranes, interfaces between compounds, or mixtures of solvents. The majority of available molecular force fields and molecular dynamics simulation software packages does not account explicitly for electronic polarization. Even the simplest charge-on-spring (COS) models have only been developed for few types of molecules. In this work, we report a polarizable COS model for cyclohexane, as this molecule is a widely used solvent, and for linear alkanes, which are also used as solvents, and are the precursors of lipids, amino acid side chains, carbohydrates, or nucleic acid backbones. The model is an extension of a nonpolarizable united-atom model for alkanes that had been calibrated against experimental values of the density, the heat of vaporization and the Gibbs free energy of hydration for each alkane. The latter quantity was used to calibrate the parameters governing the interaction of the polarizable alkanes with water. Subsequently, the model was tested for other structural, thermodynamic, dielectric, and dynamic properties such as trans/gauche ratios, excess free energy, static dielectric permittivity, and self-diffusion. A good agreement with the experimental data for a large set of properties for each considered system was obtained, resulting in a transferable set of polarizable force-field parameters for CH2, CH3, and CH4 moieties., (Copyright © 2014 Wiley Periodicals, Inc.)

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.