A set of 1630 NOE atom-atom distance upper bounds, 213 3J-coupling constant values and 200 S2 order-parameter values derived from experimental NMR data is used in molecular dynamics simulations of Hen Egg-White Lysozyme (HEWL) in aqueous solution in order to generate a Boltzmann-weighted structural ensemble for the protein in aqueous solution that is compatible with the 2043 NMR data. The two protein force fields used, the GROMOS 54A7 and 54A8 force fields, which only differ in the partial charges of the charged side chains and the chain termini of the protein, show comparable behavior. Analysis of the generated structural ensembles shows that the protein in aqueous solution adopts a greater variety of hydrogen-bond patterns than is suggested by the various X-ray crystal structures of the protein.

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.

Five sets of RDC values for the backbone of [13C,15N]-labeled Hen Egg-White Lysozyme (HEWL, 320 RDCs), obtained from NMR experiments of the protein in an ether bicelle solution at a temperature of 308 K and pH 3.8, were used to calculate RDC values by application of two methods, the alignment-tensor (AT) method and the method of magnetic-field rotational sampling (HRS), applied to five X-ray structures of HEWL, to investigate the relevance of measured RDC values for the structure determination or refinement of proteins. In contrast to other quantities Q observable by NMR, such as NOE intensities or 3J-couplings, for which a relation Q(r) between the quantity Q and a single structure r of a protein can be used to calculate average values ⟨Q(r)⟩, averaged over the Boltzmann-weighted structural ensemble of the protein at finite temperature in solution, an RDC is not defined in terms of a single structure but as an average over a slightly nonuniform rotational and orientation distribution of the protein. This averaging between large positive and negative values reduces the kHz size of a dipolar coupling (DC) by a factor of 103 to 104 to the Hz range of a residual dipolar coupling (RDC). Since the nonuniform orientation distribution can neither be measured nor faithfully mimicked at atomic resolution on a computer, RDC values for a given protein structure are commonly calculated by minimizing the difference between calculated and measured RDC values for a given set of measured target RDC values by varying the orientation distribution of the protein in one way or the other. These three features of RDCs, a very large reduction of size as a result of averaging over orientations, their definition in terms of an unknown, immeasurable orientation distribution, and their calculation using a set of target RDC values, lead to a sensitivity of the calculated RDC values to the size and type of the particular set of RDCs used in the calculation. This reduces the usefulness of measured RDCs for structure determination or refinement of proteins compared to NOE intensities or 3J-couplings.

The experimental determination of residual dipolar couplings (RDCs) rests on sampling the rotational motion of a molecule in an environment that induces a slightly nonuniform, unfortunately immeasurable, orientation distribution of the molecule in solution. Averaging over this slightly nonuniform, anisotropic distribution reduces the size of the dipolar couplings (DCs) from the kHz range to the Hz range for the resulting RDCs by a factor of 103 to 104. These features hamper the use of measured RDCs to contribute to the structure determination or refinement of (bio)molecules. The commonly used alignment-tensor (AT) methodology assumes that the immeasurable, unknown orientation distribution of the molecule can be expressed in terms of five spherical harmonic functions of order 2. Staying close to experiment, RDCs can, alternatively, be calculated from a molecular simulation by sampling the rotational motion of the molecule (MRS method) or, instead, of a vector (mfv) representing the magnetic field (HRS method). The AT and HRS methods were applied to a β-heptapeptide solvated in methanol, for which 131 NOE atom-atom distance upper bounds and 21 3J-couplings derived from NMR experiments are available and, in addition, 39 RDC values obtained for the molecule solvated in methanol with polyvinyl acetate added. In methanol at room temperature and pressure, the molecule adopts a relatively stable helical fold. It appears that MD simulation of the molecule in methanol using the GROMOS biomolecular force field already satisfies virtually all experimental data. Application of RDC restraining shows the limitations caused by the assumptions on which the AT and HRS methods rest and suggests that experimentally measured RDCs are less useful for molecular structure determination or refinement than other observable quantities that can be measured by NMR techniques. The results illustrate that in structure determination or refinement of a (bio)molecule based on experimentally measured data, it is mandatory (i) to refrain from the vacuum boundary condition and (ii) from torsional-angle restraints that do not account for the multiplicity of the inverse function of the Karplus relation expressing 3J-couplings in terms of molecular torsional angles, (iii) to allow for Boltzmann-weighted time- or molecule-averaging and, not the least, (iv) to use a force field that has an adequate basis in thermodynamic data of biomolecules.

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.

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.

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 Nat), there must be a sufficiently large number NRDC of measured RDC values available to allow the restraining forces to bias the overall orientation distribution of the molecule. If the ratio NRDC/Nat 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.

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.

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.

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.

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.

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.

This tutorial describes the practical use of some recent methodological advances implemented in the GROMOS software for biomolecular simulations. It is envisioned as a living document, with additional tutorials being added in the course of time. Currently, it consists of three distinct tutorials. The first tutorial describes the use of time-averaged restraints to enforce agreement with order parameters derived from NMR experiments. The second tutorial describes the use of extended thermodynamic integration in the double-decoupling method to compute the affinity of a small molecule to a protein. The molecule involved bears a negative charge, necessitating the application of post-simulation corrections. The third tutorial is based on the same molecular system, but computes the binding free energy from a path-sampling method with distance-field distance restraints and Hamiltonian replica exchange simulations. The tutorials are written for users with some experience in the application of molecular dynamics simulations.

Values of S2CH and S2NH 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 S2CH and S2NH order parameters towards experimentally derived target values S2CH (exp) and S2NH (exp). Application of S2CH and S2NH order-parameter restraining MD simulation to bond vectors in 63 side chains of the protein hen egg white lysozyme using 51 S2CH (exp) target values and 28 S2NH (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 S2CH order-parameter restraining of C-H bonds in methyl groups is less reliable than S2NH order-parameter restraining because of the possibly less valid assumptions and approximations used to derive experimental S2CH (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.

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.

This year, the Centre Europeen de Calcul Atomique et Moleculaire (CECAM) celebrates its 50-th anniversary. Founded in 1969 in Orsay near Paris, it later moved to Lyon and in 2008 to Lausanne. It is an organization devoted to the promotion of fundamental research on advanced computational methods and their application in condensed matter science. Its main vehicle to this end is the organization of workshops. The key role of an eight-week workshop held forty-three years ago, characterized by an open exchange of scientific ideas and a foresight regarding the topics relevant to a proper dynamic simulation of bio-molecules such as proteins, is remembered, together with the issues discussed at the time. These are still relevant today.

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 Qexp of experimentally observable quantities Q. One must consider 1) the accuracy of Qexp , 2) the accuracy of the function Q(rN ) used to calculate a Q-value based on a molecular configuration rN of N particles, 3) the sensitivity of the function Q(rN ) to the configuration rN , 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.

The use of a supra-molecular coarse-grained (CG) model for liquid water as solvent in molecular dynamics simulations of biomolecules represented at the fine-grained (FG) atomic level of modelling may reduce the computational effort by one or two orders of magnitude. However, even if the pure FG model and the pure CG model represent the properties of the particular substance of interest rather well, their application in a hybrid FG/CG system containing varying ratios of FG versus CG particles is highly non-trivial, because it requires an appropriate balance between FG-FG, FG-CG, and CG-CG energies, and FG and CG entropies. Here, the properties of liquid water are used to calibrate the FG-CG interactions for the simple-point-charge water model at the FG level and a recently proposed supra-molecular water model at the CG level that represents five water molecules by one CG bead containing two interaction sites. Only two parameters are needed to reproduce different thermodynamic and dielectric properties of liquid water at physiological temperature and pressure for various mole fractions of CG water in FG water. The parametrisation strategy for the FG-CG interactions is simple and can be easily transferred to interactions between atomistic biomolecules and CG water.