A Method to Derive Structural Information on Molecules from Residual Dipolar Coupling NMR 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 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.

On the use of intra-molecular distance and angle constraints to lengthen the time step in molecular and stochastic dynamics simulations of proteins

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.

On the use of multiple-time-step algorithms to save computing effort in molecular dynamics simulations of proteins

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.

A method to apply bond-angle constraints in molecular dynamics simulations

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.

On the use of 3J-coupling NMR data to derive structural information on proteins

Values of ³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 ³J-coupling values (4-8 Hz) to the averaging period that must cover the conformational variability of the torsional angle related to the ³J-coupling, and due to the difficulty of handling the multiple-valued character of the inverse Karplus relation between torsional angle and ³J-coupling. Both problems can be solved by using ³J-coupling time-averaging local-elevation restraining MD simulation. Application to the protein hen egg white lysozyme using 213 backbone and side-chain ³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 ³J-coupling data in protein structure determination.

On the Effect of the Various Assumptions and Approximations used in Molecular Simulations on the Properties of Bio-Molecular Systems: Overview and Perspective on Issues

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.

On the Use of Side-Chain NMR Relaxation Data to Derive Structural and Dynamical Information on Proteins: A Case Study Using Hen Lysozyme

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.

Algorithms to apply dihedral-angle constraints in molecular or stochastic dynamics simulations

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.

Conformational Properties of the Chemotherapeutic Drug Analogue Epothilone A: How to Model a Flexible Protein Ligand Using Scarcely Available Experimental Data

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.

Validation of Molecular Simulation: An Overview of Issues

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.

Using Complementary NMR Data Sets To Detect Inconsistencies and Model Flaws in the Structure Determination of Human Interleukin-4

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² order parameters, and ³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.

Interpretation of Seemingly Contradictory Data: Low NMR S2 Order Parameters Observed in Helices and High NMR S2 Order Parameters in Disordered Loops of the Protein hGH at Low pH

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² 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² order parameter restraining MD simulations, this paradox has been resolved. Low S² values in helices are due to the presence of a mixture of 3₁₀ -helical and α-helical hydrogen bonds. High S² 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.

Deriving Structural Information from Experimentally Measured Data on Biomolecules

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.

On the use of time-averaging restraints when deriving biomolecular structure from ³J -coupling values obtained from NMR experiments

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.

Investigation of the structural preference and flexibility of the loop residues in amyloid fibrils of the HET-s prion

The structural variability of HET-s(218–289) loops is restricted by the β-sheet core.

Rapid Sampling of Folding Equilibria of β-Peptides in Methanol Using a Supramolecular Solvent Model

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.

Computational Analysis of the Mechanism and Thermodynamics of Inhibition of Phosphodiesterase 5A by Synthetic Ligands

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

Theoretical Investigation of Solvent Effects on Glycosylation Reactions: Stereoselectivity Controlled by Preferential Conformations of the Intermediate Oxacarbenium-Counterion Complex

The mechanism of solvent effects on the stereoselectivity of glycosylation reactions is investigated using quantum-mechanical (QM) calculations and molecular dynamics (MD) simulations, considering a methyl-protected glucopyranoside triflate as a glycosyl donor equivalent and the solvents acetonitrile, ether, dioxane, or toluene, as well as gas-phase conditions (vacuum). The QM calculations on oxacarbenium-solvent complexes do not provide support to the usual solvent-coordination hypothesis, suggesting that an experimentally observed β-selectivity (α-selectivity) is caused by the preferential coordination of a solvent molecule to the reactive cation on the α-side (β-side) of the anomeric carbon. Instead, explicit-solvent MD simulations of the oxacarbenium-counterion (triflate ion) complex (along with corresponding QM calculations) are compatible with an alternative mechanism, termed here the conformer and counterion distribution hypothesis. This new hypothesis suggests that the stereoselectivity is dictated by two interrelated conformational properties of the reactive complex, namely, (1) the conformational preferences of the oxacarbenium pyranose ring, modulating the steric crowding and exposure of the anomeric carbon toward the α or β face, and (2) the preferential coordination of the counterion to the oxacarbenium cation on one side of the anomeric carbon, hindering a nucleophilic attack from this side. For example, in acetonitrile, the calculations suggest a dominant B2,5 ring conformation of the cation with preferential coordination of the counterion on the α side, both factors leading to the experimentally observed β selectivity. Conversely, in dioxane, they suggest a dominant (4)H3 ring conformation with preferential counterion coordination on the β side, both factors leading to the experimentally observed α selectivity.

Flexible Boundaries for Multiresolution Solvation: An Algorithm for Spatial Multiscaling in Molecular Dynamics Simulations

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.

On the use of a weak-coupling thermostat in replica-exchange molecular dynamics 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.

An improved simple polarisable water model for use in biomolecular simulation

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.

Supra-Atomic Coarse-Grained GROMOS Force Field for Aliphatic Hydrocarbons in the Liquid Phase

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.

A comparison of pathway-independent and pathway-dependent methods in the calculation of conformational free enthalpy differences

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.

Time-averaged order parameter restraints in molecular dynamics simulations

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.

Practical Aspects of Free-Energy Calculations: A Review

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.

A polarizable empirical force field for molecular dynamics simulation of liquid hydrocarbons

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.

Exploration of swapping enzymatic function between two proteins: a simulation study of chorismate mutase and isochorismate pyruvate lyase

The enzyme chorismate mutase EcCM from Escherichia coli catalyzes one of the few pericyclic reactions in biology, the transformation of chorismate to prephenate. The isochorismate pyruvate lyase PchB from Pseudomonas aeroginosa catalyzes another pericyclic reaction, the isochorismate to salicylate transformation. Interestingly, PchB possesses weak chorismate mutase activity as well thus being able to catalyze two distinct pericyclic reactions in a single active site. EcCM and PchB possess very similar folds, despite their low sequence identity. Using molecular dynamics simulations of four combinations of the two enzymes (EcCM and PchB) with the two substrates (chorismate and isochorismate) we show that the electrostatic field due to EcCM at atoms of chorismate favors the chorismate to prephenate transition and that, analogously, the electrostatic field due to PchB at atoms of isochorismate favors the isochorismate to salicylate transition. The largest differences between EcCM and PchB in electrostatic field strengths at atoms of the substrates are found to be due to residue side chains at distances between 0.6 and 0.8 nm from particular substrate atoms. Both enzymes tend to bring their non-native substrate in the same conformation as their native substrate. EcCM and to a lower extent PchB fail in influencing the forces on and conformations of the substrate such as to favor the other chemical reaction (isochorismate pyruvate lyase activity for EcCM and chorismate mutase activity for PchB). These observations might explain the difficulty of engineering isochorismate pyruvate lyase activity in EcCM by solely mutating active site residues.

On the behavior of water at subfreezing temperatures in a protein crystal: evidence of higher mobility than in bulk water

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.

Influence of variation of a side chain on the folding equilibrium of a β-peptide: limitations of one-step perturbation

In a recent study (Lin et al., Helv. Chim. Acta 2011, 94, 597), the one-step perturbation method was applied to tackle a challenging computational problem, that is, the calculation of the folding free enthalpies ΔGF,U of six hepta-β-peptides with different, Ala, Val, Leu, Ile, Ser, or Thr, side chains in the fifth residue. The ΔGF,U values obtained using one-step perturbation based on a single molecular dynamics simulation of a judiciously chosen reference state with soft-core atoms in the side chain of the fifth residue showed an overall accuracy of about kB T for the four peptides with nonpolar side chains, but twice as large deviations were observed for the peptides with polar side chains. Here, alternative reference-state Hamiltonians that better cover the conformational space relevant to these peptides are investigated, and post simulation rotational sampling of the χ1 and χ2 torsional angles of the fifth residue is carried out to sample different orientations of the side chain. A reference state with rather soft atoms yields accurate ΔGF,U values for the peptides with the Ser and Thr side chains, but it failed to correctly predict the folding free enthalpy for one peptide with a nonpolar side chain, that is, Leu. Based on the results and those of earlier studies, possible ways to improve the accuracy of the efficient one-step perturbation technique to compute free enthalpies of folding are discussed.

Molecular dynamics simulations of barley and maize lipid transfer proteins show different ligand binding preferences in agreement with experimental data

Experimental studies of barley and maize lipid transfer proteins (LTPs) show that the two proteins bind the ligand palmitate in opposite orientations in their internal cavities. Moreover, maize LTP is reported to bind the ligand caprate in the internal cavity in a mixture of two orientations with approximately equal occupancy. Six 30 ns molecular dynamics (MD) simulations of maize and barley LTP with ligands bound in two orientations (modes M and B) have been used to understand the different ligand binding preferences. The simulations show that both maize and barley LTP could bind palmitate in the orientation observed experimentally for maize LTP (mode M), with the predominant interaction being a salt bridge between the ligand carboxylate headgroup and a conserved arginine side chain. However, the simulation of barley LTP with palmitate in the mode B orientation shows the most favorable protein-ligand interaction energy. In contrast, the simulations of maize LTP with palmitate and with caprate in the mode B orientation show no persistent ligand binding, the ligands leaving the cavity during the simulations. Sequence differences between maize and barley LTP in the AB loop region, in residues at the base of the hydrophobic cavity, and in the helix A region are identified as contributing to the different behavior. The simulations reproduce well the experimentally observed binding preferences for palmitate and suggest that the experimental data for maize LTP with caprate reflect ligand mobility in binding mode M rather than the population of binding modes M and B.