Search

The level of accuracy that can be achieved by a force field is influenced by choices made in the interaction-function representation and in the relevant simulation parameters. These choices, referred to here as functional-form variants (FFVs), include for example the model resolution, the charge-derivation procedure, the van der Waals combination rules, the cutoff distance, and the treatment of the long-range interactions. Ideally, assessing the effect of a given FFV on the intrinsic accuracy of the force-field representation requires that only the specific FFV is changed and that this change is performed at an optimal level of parametrization, a requirement that may prove extremely challenging to achieve in practice. Here, we present a first attempt at such a comparison for one specific FFV, namely the choice of a united-atom (UA) versus an all-atom (AA) resolution in a force field for saturated acyclic (halo)alkanes. Two force-field versions (UA vs AA) are optimized in an automated way using the CombiFF approach against 961 experimental values for the pure-liquid densities ρliq and vaporization enthalpies ΔHvap of 591 compounds. For the AA force field, the torsional and third-neighbor Lennard–Jones parameters are also refined based on quantum-mechanical rotational-energy profiles. The comparison between the UA and AA resolutions is also extended to properties that have not been included as parameterization targets, namely the surface-tension coefficient γ, the isothermal compressibility κT, the isobaric thermal-expansion coefficient αP, the isobaric heat capacity cP, the static relative dielectric permittivity ϵ, the self-diffusion coefficient D, the shear viscosity η, the hydration free energy ΔGwat, and the free energy of solvation ΔGche in cyclohexane. For the target properties ρliq and ΔHvap, the UA and AA resolutions reach very similar levels of accuracy after optimization. For the nine other properties, the AA representation leads to more accurate results in terms of η; comparably accurate results in terms of γ, κT, αP, ϵ, D, and ΔGche; and less accurate results in terms of cP and ΔGwat. This work also represents a first step toward the calibration of a GROMOS-compatible force field at the AA resolution., Journal of Chemical Theory and Computation, 18 (11), ISSN:1549-9618, ISSN:1549-9626

The calibration of torsional interaction terms by fitting relative gas-phase conformational energies against their quantum-mechanical values is a common procedure in force-field development. However, much less attention has been paid to the optimization of third-neighbor nonbonded interaction parameters, despite their strong coupling with the torsions. This article introduces an algorithm termed LLS-SC, aimed at simultaneously parametrizing torsional and third-neighbor interaction terms based on relative conformational energies. It relies on a self-consistent (SC) procedure where each iteration involves a linear least-squares (LLS) regression followed by a geometry optimization of the reference structures. As a proof-of-principle, this method is applied to obtain torsional and third-neighbor interaction parameters for aliphatic chains in the context of the GROMOS 53A6 united-atom force field. The optimized parameter set is compared to the original one, which has been fitted manually against thermodynamic properties for small linear alkanes. The LLS-SC implementation is freely available under http://github.com/mssm-labmmol/profiler. ISSN:0192-8651 ISSN:1096-987X

Alchemical free energy calculations generally require intermediate states along a coupling parameter lambda to establish sufficient phase space overlap for obtaining converged results. Such intermediate states can also be engineered to lower the energy barriers and, consequently, reduce the required sampling time. The recently introduced lambda-enveloping distribution sampling (lambda-EDS) scheme combines the properties of the minimum variance pathway and the EDS methods to improve sampling and allow for larger steps along the alchemical pathway compared to conventional approaches. This scheme also eliminates the need for soft-core potentials and retains the behavior of conventional lambda-intermediate states as a limiting case. In this study, an automated procedure is developed to select the parameters of lambda-EDS for optimal performance. The underlying theory is illustrated based on simulations of simple test systems (bond length changes in harmonic oscillators, mutations of dihedral angles, and charge creation in water), as well as on the calculation of the absolute hydration free energies of 12 small organic molecules. © 2021 American Chemical Society ISSN:1549-9618 ISSN:1549-9626

Free-energy differences between pairs of end-states can be estimated based on molecular dynamics (MD) simulations using standard pathway-dependent methods such as thermodynamic integration (TI), free-energy perturbation, or Bennett's acceptance ratio. Replica-exchange enveloping distribution sampling (RE-EDS), on the other hand, allows for the sampling of multiple end-states in a single simulation without the specification of any pathways. In this work, we use the RE-EDS method as implemented in GROMOS together with generalized AMBER force-field (GAFF) topologies, converted to a GROMOS-compatible format with a newly developed GROMOS++ program amber2gromos, to compute relative hydration free energies for a series of benzene derivatives. The results obtained with RE-EDS are compared to the experimental data as well as calculated values from the literature. In addition, the estimated free-energy differences in water and in vacuum are compared to values from TI calculations carried out with GROMACS. The hydration free energies obtained using RE-EDS for multiple molecules are found to be in good agreement with both the experimental data and the results calculated using other free-energy methods. While all considered free-energy methods delivered accurate results, the RE-EDS calculations required the least amount of total simulation time. This work serves as a validation for the use of GAFF topologies with the GROMOS simulation package and the RE-EDS approach. Furthermore, the performance of RE-EDS for a large set of 28 end-states is assessed with promising results., Journal of Chemical Information and Modeling, 62 (12), ISSN:1549-9596, ISSN:0095-2338, ISSN:1520-5142

Molecular dynamics studies have demonstrated that molecular water at an interface, with either a gas or a solid, displays anisotropic orientational behavior in contrast to its bulk counterpart. This effect has been recently implicated in the like-charge attraction problem for colloidal particles in solution. Here, negatively charged particles in solution display a long-ranged attraction where continuum electrostatic theory predicts monotonically repulsive interactions, particularly in solutions with monovalent salt ions at low ionic strength. Anisotropic orientational behavior of solvent molecules at an interface gives rise to an excess interfacial electrical potential which we suggest generates an additional solvation contribution to the total free energy that is traditionally overlooked in continuum descriptions of interparticle interactions in solution. In the present investigation we perform molecular dynamics simulation based calculations of the interfacial potential using realistic surface models representing various chemistries as well as different solvents. Similar to previous work that focused on simple model surfaces constructed by using oxygen atoms, we find that solvents at more realistic model surfaces exhibit substantial anisotropic orientational behavior. We explore the dependence of the interfacial solvation potential on surface properties such as surface group chemistry and group density at silica and carboxylated polystyrene interfaces. For water, we note surprisingly good agreement between results obtained for a simple O-atom wall and more complex surface models, suggesting a general qualitative consistency of the interfacial solvation effect for surfaces in contact with water. In contrast, for an aprotic solvent such as DMSO, surface chemistry appears to exert a stronger influence on the sign and magnitude of the interfacial solvation potential. The study carries broad implications for molecular-scale interactions and may find relevance in explaining a range of phenomena in soft-matter physics and cell biology., The Journal of Physical Chemistry B, 126 (25), ISSN:1520-6106, ISSN:1520-5207, ISSN:1089-5647

The CombiFF approach is a workflow for the automated refinement of force-field parameters against experimental condensed-phase data, considering entire classes of organic molecules constructed using a fragment library via combinatorial isomer enumeration. One peculiarity of this approach is that it relies on an electronegativity-equalization scheme to account for induction effects within molecules, with values of the atomic hardness and electronegativity as electrostatic parameters, rather than the partial charges themselves. In a previous article [M. P. Oliveira, M. Andrey, S. R. Rieder, L. Kern, D. F. Hahn, S. Riniker, B. A. C. Horta and P. H. Hünenberger, J. Chem. Theory. Comput. 2020, 16, 7525], CombiFF was introduced and applied to calibrate a GROMOS-compatible united-atom force field for the saturated acyclic (halo-)alkane family. Here, this scheme is employed for the construction of a corresponding force field for saturated acyclic compounds encompassing eight common chemical functional groups involving oxygen and/or nitrogen atoms, namely: ether, aldehyde, ketone, ester, alcohol, carboxylic acid, amine, and amide. Monofunctional as well as homo-polyfunctional compounds are considered. A total of 1712 experimental liquid densities ρliq and vaporization enthalpies ΔHvap concerning 1175 molecules are used for the calibration (339 molecules) and validation (836 molecules) of the 102 non-bonded interaction parameters of the force field. Using initial parameter values based on the GROMOS 2016H66 parameter set, convergence is reached after five iterations. Given access to one processor per simulated system, this operation only requires a few days of wall-clock computing time. After optimization, the root-mean-square deviations from experiment are 29.9 (22.4) kg m−3 for ρliq and 4.1 (5.5) kJ mol−1 for ΔHvap for the calibration (validation) set. Thus, a very good level of agreement with experiment is achieved in terms of these two properties, although the errors are inhomogeneously distributed across the different chemical functional groups., Physical Chemistry Chemical Physics, 23 (33), ISSN:1463-9084, ISSN:1463-9076

A new approach termed Adaptive Solvent-Scaling (AdSoS) is introduced for performing simulations of a solute embedded in a fine-grained (FG) solvent region itself surrounded by a coarse-grained (CG) solvent region, with a continuous FG ↔ CG switching of the solvent resolution across a buffer layer. Instead of relying on a distinct CG solvent model, the AdSoS scheme is based on CG models defined by a dimensional scaling of the FG solvent by a factor s, accompanied by an s-dependent modulation of the atomic masses and interaction parameters. The latter changes are designed to achieve an isomorphism between the dynamics of the FG and CG models, and to preserve the dispersive and dielectric solvation properties of the solvent with respect to a solute at FG resolution. This scaling approach offers a number of advantages compared to traditional coarse-graining: (i) the CG parameters are immediately related to those of the FG model (no need to parameterize a distinct CG model); (ii) nearly ideal mixing is expected for CG variants with similar s-values (ideal mixing holding in the limit of identical s-values); (iii) the solvent relaxation timescales should be preserved (no dynamical acceleration typical for coarse-graining); (iv) the graining level N

Replica-exchange enveloping distribution sampling (RE-EDS) is a pathway-independent multistate free-energy method currently implemented in the GROMOS software package for molecular dynamics (MD) simulations. It has a high intrinsic sampling efficiency as the interactions between the unperturbed particles have to be calculated only once for multiple end-states. As a result, RE-EDS is an attractive method for the calculation of relative solvation and binding free energies. An essential requirement for reaching this high efficiency is the separability of the nonbonded interactions into solute–solute, solute–environment, and environment–environment contributions. Such a partitioning is trivial when using a Coulomb term with a reaction-field (RF) correction to model the electrostatic interactions but not when using lattice-sum schemes. To avoid cutoff artifacts, the RF correction is typically used in combination with a charge-group-based cutoff, which is not supported by most small-molecule force fields as well as other MD engines. To address this issue, we investigate the combination of RE-EDS simulations with a recently introduced RF scheme including a shifting function that enables the rigorous calculation of RF electrostatics with atom-based cutoffs. The resulting approach is validated by calculating solvation free energies with the generalized AMBER force field in water and chloroform using both the GROMOS software package and a proof-of-concept implementation in OpenMM., The Journal of Chemical Physics, 157 (10), ISSN:0021-9606, ISSN:1089-7690

The effect of different treatments of the nonbonded interactions in simulations employing the recently introduced GROMOS-compatible 2016H66 force field is evaluated based on calculations carried out with the GROMACS software. This is done considering four thermodynamic and transport properties (pure liquid density, vaporization enthalpy, surface-tension coefficient, and self-diffusion constant) of 58 organic liquids representative of the chemical groups alcohol, ether, aldehyde, ketone, carboxylic acid, ester, amine, amide, thiol, sulfide, disulfide, and aromatic compounds, also including water (SPC model). A dipalmitoylphosphatidylcholine bilayer system is considered as well. The simulated properties are found to be very sensitive to the treatment of the long-range dispersion interactions, notably for the least polar systems. In general, the treatment of the long-range electrostatic or Lennard-Jones interactions using homogeneous correction terms or lattice-sum approaches yield similar results, with punctual discrepancies. The combination of a lattice-sum approach for the electrostatic interactions with a straight-cutoff truncation of the Lennard-Jones interactions at a distance of at least 1.2 nm is found to represent a good compromise setup within GROMACS for achieving compatibility with the reference results obtained using GROMOS as well as a comparable level of agreement with the experimental data. This study also reveals two potential issues with the GROMACS software, related to an incorrect calculation of the pressure when using LINCS in version 4.0.7 and an inadequate implementation of the twin-range scheme in version 5.1.2.

Experimental solvation free energies are nowadays commonly included as target properties in the validation of condensed-phase force fields, sometimes even in their calibration. In a previous article [Kashefolghetaet al.,J. Chem. Theory. Comput., 2020,16, 7556-7580], we showed how the involved comparison between experimental and simulation results could be made more systematic by considering a full matrix of cross-solvation free energies. For a set ofNmolecules that are all in the liquid state under ambient conditions, such a matrix encompassesN×Nentries for considering each of theNmolecules either as solute (A) or as solvent (B). In the quoted study, a cross-solvation matrix of 25 × 25 experimental value was introduced, considering 25 small molecules representative for alkanes, chloroalkanes, ethers, ketones, esters, alcohols, amines, and amides. This experimental data was used to compare the relative accuracies of four popular condensed-phase force fields, namely GROMOS-2016H66, OPLS-AA, AMBER-GAFF, and CHARMM-CGenFF. In the present work, the comparison is extended to five additional force fields, namely GROMOS-54A7, GROMOS-ATB, OPLS-LBCC, AMBER-GAFF2, and OpenFF. Considering these nine force fields, the correlation coefficients between experimental values and simulation results range from 0.76 to 0.88, the root-mean-square errors (RMSEs) from 2.9 to 4.8 kJ mol−1, and average errors (AVEEs) from −1.5 to +1.0 kJ mol−1. In terms of RMSEs, GROMOS-2016H66 and OPLS-AA present the best accuracy (2.9 kJ mol−1), followed by OPLS-LBCC, AMBER-GAFF2, AMBER-GAFF, and OpenFF (3.3 to 3.6 kJ mol−1), and then by GROMOS-54A7, CHARM-CGenFF, and GROMOS-ATB (4.0 to 4.8 kJ mol−1). These differences are statistically significant but not very pronounced, and are distributed rather heterogeneously over the set of compounds within the different force fields., Physical Chemistry Chemical Physics, 23 (23), ISSN:1463-9084, ISSN:1463-9076

Direct optimization against experimental condensed-phase properties concerning small organic molecules still represents the most reliable way to calibrate the empirical parameters of a force field. However, compared to a corresponding calibration against quantum-mechanical (QM) calculations concerning isolated molecules, this approach is typically very tedious and time-consuming. The present article describes an integrated scheme for the automated refinement of force-field parameters against experimental condensed-phase data, considering entire classes of organic molecules constructed using a fragment library via combinatorial isomer enumeration. The main steps of the scheme, referred to as CombiFF, are as follows: (i) definition of a molecule family; (ii) combinatorial enumeration of all isomers; (iii) query for experimental data; (iv) automatic construction of the molecular topologies by fragment assembly; and (v) iterative refinement of the force-field parameters considering the entire family. As a first application, CombiFF is used here to design a GROMOS-compatible united-atom force field for the saturated acyclic haloalkane family. This force field relies on an electronegativity-equalization scheme for the atomic partial charges and involves no specific terms for σ-holes and halogen bonding. A total of 749 experimental liquid densities ρliq and vaporization enthalpies ΔHvap concerning 486 haloalkanes are considered for calibration and validation. The resulting root-mean-square deviations from experiment are 49.8 (27.6) kg·m–3 for ρliq and 2.7 (1.8) kJ·mol–1 for ΔHvap for the calibration (validation) set. The values are lower for the validation set which contains larger molecules (stronger influence of purely aliphatic interactions). The trends in the optimized parameters along the halogen series and across the compound family are in line with chemical intuition based on considerations related to size, polarizability, softness, electronegativity, induction, and hyperconjugation. This observation is particularly remarkable considering that the force-field calibration did not involve any QM calculation. Once the time-consuming task of target-data selection/curation has been performed, the optimization of a force field only takes a few days. As a result, CombiFF enables an easy assessment of the consequences of functional-form decisions on the accuracy of a force field at an optimal level of parametrization., Journal of Chemical Theory and Computation, 16 (12), ISSN:1549-9618, ISSN:1549-9626

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.

Alchemical free energy calculations typically rely on intermediate states to bridge between the relevant phase spaces of the two end states. These intermediate states are usually created by mixing the energies or parameters of the end states according to a coupling parameter λ. The choice of the procedure has a strong impact on the efficiency of the calculation, as it affects both the encountered energy barriers and the phase space overlap between the states. The present work builds on the connection between the minimum variance pathway (MVP) and enveloping distribution sampling (EDS). It is shown that both methods can be regarded as special cases of a common scheme referred to as λ-EDS, which can also reproduce the behavior of conventional λ-intermediate states. A particularly attractive feature of λ-EDS is its ability to emulate the use of soft core potentials (SCP) while avoiding the associated computational overhead when applying efficient free energy estimators such as the multistate Bennett’s acceptance ratio (MBAR). The method is illustrated for both relative and absolute free energy calculations considering five benchmark systems. The first two systems (charge inversion and cavity creation in a dipolar solvent) demonstrate the use of λ-EDS as an alternative coupling scheme in the context of thermodynamic integration (TI). The three other systems (change of bond length, change of dihedral angles, and cavity creation in water) investigate the efficiency and optimal choice of parameters in the context of free energy perturbation (FEP) and Bennett’s acceptance ratio (BAR). It is shown that λ-EDS allows larger steps along the alchemical pathway than conventional intermediate states. © 2020 American Chemical Society. ISSN:1549-9596 ISSN:0095-2338 ISSN:1520-5142

In molecular dynamics (MD) simulations of condensed-phase systems, straight-cutoff truncation of the non-bonded interactions is well known to cause cutoff noise and serious artifacts in many simulated properties. These effects can be drastically reduced by applying the truncation based on distances between neutral charge groups (CG) rather than between individual atoms (AT). In addition, the mean effect of the omitted electrostatic interactions beyond the cutoff distance can be reintroduced using the reaction-field (RF) method, where the medium outside the cutoff sphere is approximated as a dielectric continuum of permittivity equal to that of the solvent. The RF scheme is generally applied with CG truncation. This is justified for low solvent permittivities, where the RF correction is small and an AT truncation would lead to severe issues, just as in the straight-cutoff case. However, it is less appropriate for solvents with high permittivities, where the RF correction acts as a physically motivated shifting function, and a CG truncation may in turn lead to artifacts and poorer energy conservation. In this study, we assess the impact of truncation artifacts considering the 57 organic liquids which were used in the calibration of the GROMOS-compatible 2016H66 force field. Combinations of shifting or switching schemes with RF-based electrostatic interactions as well as van der Waals (Lennard-Jones) interactions are then introduced to resolve the issues with AT truncation. These shifting and switching schemes have the following properties: (i) they bring the force but not the potential energy to zero at the cutoff; (ii) as a result, they lead to a modification of the interaction that is comparatively small; (iii) they permit to conduct rigorously conservative simulations; (iv) the energies can easily be corrected back to the unmodified form, either on the fly or in a post-processing step. The mathematical formalism of these schemes is presented in detail, and their validation is performed using the 57 organic liquids., Physical Chemistry Chemical Physics, 22 (45), ISSN:1463-9084, ISSN:1463-9076

Experimental solvation free energies are nowadays commonly included as target properties in the validation and sometimes even in the calibration of condensed-phase force fields. However, this is often done in a nonsystematic fashion, by considering available solvation free energies involving an arbitrary collection of solutes in a limited set of solvents (e.g., water, octanol, chloroform, cyclohexane, or hexane). Here, this approach is made more systematic by introducing the concept of cross-solvation free energies ΔsGA:B⊖ for a set of N molecules that are all in the liquid state under ambient conditions, namely the matrix of N2 entries for ΔsGA:B⊖ considering each of the N molecules either as a solute (A) or as a solvent (B). Relying on available experimental literature followed by careful data curation, a complete ΔsGA:B⊖ matrix of 625 entries is constructed for 25 molecules with one to seven carbon atoms representative for alkanes, chloroalkanes, ethers, ketones, esters, alcohols, amines, and amides. This matrix is then used to compare the relative accuracies of four popular condensed-phase force fields: GROMOS-2016H66, OPLS-AA, AMBER-GAFF, and CHARMM-CGenFF. In broad terms, and in spite of very different force-field functional-form choices and parametrization strategies, the four force fields are found to perform similarly well. Relative to the experimental values, the root-mean-square errors range between 2.9 and 4.0 kJ·mol–1 (lowest value of 2.9 for GROMOS and OPLS), and the average errors range between −0.8 and +1.0 kJ·mol–1 (lowest magnitude of 0.2 for AMBER and CHARMM). These differences are statistically significant but not very pronounced, especially considering the influence of outliers, some of which possibly caused by inaccurate experimental data., Journal of Chemical Theory and Computation, 16 (12), ISSN:1549-9618, ISSN:1549-9626

We recently introduced a method called conveyor belt (CB) thermodynamic integration (TI) for the calculation of alchemical free-energy differences based on molecular dynamics simulations. In the present work, the CBTI approach is generalized to conformational free-energy changes, i.e., to the determination of the potential of mean force (PMF) along a conformational coordinate ξ of interest. The proposed conveyor belt umbrella sampling (CBUS) scheme relies on the parallel simulation of K replicas k = 0,1, ..., K - 1 of the system, with K even. For each replica k, the instantaneous value of ξ is restrained to an anchor value λk. The latter anchor points are equally spaced along a forward-turn-backward-turn path (i.e., a CB) between two extreme values defining the ξ-range of interest. The rotation of the CB is controlled by a variable Λ (range from 0 to 2π) which evolves dynamically along the simulation. The evolution of Λ results from the forces exerted by the restraining potentials on the anchor points, taken equal and opposite to those they exert on the replicas. Because these forces tend to cancel out along the CB, the dynamics of Λ is essentially diffusive, and the continuous distribution of ξ-values sampled by the replica system is automatically close to homogeneous. The latter feature represents an advantage over direct counting (DCNT) and traditional umbrella sampling (TRUS), shared to some extent with replica-exchange umbrella sampling (REUS). In this work, the CBUS scheme is introduced and compared to the three latter schemes in the calculation of 45 standard absolute binding free energies. These correspond to the binding of five alkali cations to three crown ethers in three solvents. Different free-energy estimators are considered for the PMF calculation, and the calculated values are also compared to those of a previous study relying on an alchemical path, as well as to experimental data.

Alchemical free energy calculations using conventional molecular dynamics and thermodynamic integration rely on simulations performed at fixed values of the coupling parameter λ. When multiple conformers in equilibrium are separated by high barriers in the space orthogonal to λ, proper convergence may require extremely long simulations. Four main strategies can be employed to address this orthogonal-sampling problem: (a) λ-variations, where λ can change along the simulations to circumvent barriers; (b) λ-extrapolations, where statistical information is transferred between λ-points; (c) specific biasing, where orthogonal barriers are reduced using a biasing potential designed specifically for the system; and (d) generic biasing, where orthogonal barriers are reduced using a generic approach. Here, we investigate the relative merits of the first three strategies considering two benchmark systems. The KXK system involves a mutation of the central residue in a tripeptide to a glycine and the XTP system involves a hydrogen-to-bromine mutation in the base of a nucleotide. Three sampling methods are compared, the latter two involving λ-variations: molecular dynamics simulations with fixed λ-points, Hamiltonian replica exchange, and the recently introduced conveyor belt method. Two free energy estimators are applied, the second one involving λ-extrapolations: thermodynamic integration with Simpson quadrature and the multistate Bennett acceptance ratio. Finally, three different seeding schemes are considered for the generation of the initial configurations. For both benchmark systems, λ-extrapolations are found to provide little gain, whereas λ-variations can significantly enhance the convergence. They are sufficient on their own if the orthogonal barriers are low in at least one state (e.g., the glycine state in KXK). However, if the orthogonal barriers are high over the entire λ-range (e.g., the XTP system), λ-variations are only effective when applied together with a specific biasing for introducing such a low-barrier state.

A recently proposed thermodynamic integration (TI) approach formulated in the framework of quantum mechanical/molecular mechanical molecular dynamics (QM/MM MD) simulations is applied to study the ...

Polyoxyethylene glycol alkyl ether amphiphiles (CiEj) are important nonionic surfactants, often used for biophysical and membrane protein studies. In this work, we extensively test the GROMOS-compatible 2016H66 force field in molecular dynamics simulations involving the lamellar phase of a series of CiEj surfactants, namely C12E2, C12E3, C12E4, C12E5, and C14E4. The simulations reproduce qualitatively well the monitored structural properties and their experimental trends along the surfactant series, although some discrepancies remain, in particular in terms of the area per surfactant, the equilibrium phase of C12E5, and the order parameters of C12E3, C12E4, and C12E5. The polar head of the CiEj surfactants is highly hydrated, almost like a single polyethyleneoxide (PEO) molecule at full hydration, resulting in very compact conformations. Within the bilayer, all CiEj surfactants flip-flop spontaneously within tens of nanoseconds. Water-permeation is facilitated, and the bending rigidity is 4 to 5 times low...

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

The asymmetric solvation of ions can be defined as the tendency of a solvent to preferentially solvate anions over cations or cations over anions, at identical ionic charge magnitudes and effective sizes. Taking water as a reference, these effects are quantified experimentally for many solvents by the relative acity (A) and basity (B) parameters of the Swain scale. The goal of the present study is to investigate the asymmetric solvation of ions using molecular dynamics simulations, and to connect the results to this empirical scale. To this purpose, the charging free energies of alkali and halide ions, and of their hypothetical oppositely charged counterparts, are calculated in a variety of solvents. In a first set of calculations, artificial solvent models are considered that present either a charge or a shape asymmetry at the molecular level. The solvation asymmetry, probed by the difference in charging free energy between the two oppositely charged ions, is found to encompass a term quadratic in the ion charge, related to the different solvation structures around the anion and cation, and a term linear in the ion charge, related to the solvation structure around the uncharged ion-sized cavity. For these simple solvent models, the two terms are systematically counteracting each other, and it is argued that only the quadratic term should be retained when comparing the results of simulations involving physical solvents to experimental data. In a second set of calculations, 16 physical solvents are considered. The theoretical estimates for the acity A are found to correlate very well with the Swain parameters, whereas the correlation for B is very poor. Based on this observation, the Swain scale is reformulated into a new scale involving an asymmetry parameter Σ, positive for acitic solvents and negative for basitic ones, and a polarity parameter Π. This revised scale has the same predictive power as the original scale, but it characterizes asymmetry in an absolute sense, the atomistic simulations playing the role of an extra-thermodynamic assumption, and is optimally compatible with the simulation results. Considering the 55 solvents in the Swain set, it is observed that a moderate basity (Σ between -0.9 and -0.3, related to electronic polarization) represents the baseline for most solvents, while a highly variable acity (Σ between 0.0 and 3.0, related to hydrogen-bond donor capacity modulated by inductive effects) represents a landmark of protic solvents.

This article describes a revised version 56A6(CARBO_R) of the GROMOS 56A6(CARBO) force field for hexopyranose-based carbohydrates. The simulated properties of unfunctionalized hexopyranoses are unaltered with respect to 56A6CARBO . In the context of both O1 -alkylated hexopyranoses and oligosaccharides, the revision stabilizes the regular (4) C1 chair for α-anomers, with the opposite effect for β-anomers. As a result, spurious ring inversions observed in α(1→4)-linked chains when using the original 56A6(CARBO) force field are alleviated. The (4) C1 chair is now the most stable conformation for all d-hexopyranose residues, irrespective of the linkage type and anomery, and of the position of the residue along the chain. The methylation of a d-hexopyranose leads to a systematic shift in the ring-inversion free energy ((4) C1 to (1) C4 ) by 7-8 kJ mol(-1), positive for the α-anomers and negative for the β-anomers, which is qualitatively compatible with the expected enhancement of the anomeric effect upon methylation at O1. The ring-inversion free energies for residues within chains are typically smaller in magnitude compared to those of the monomers, and correlate rather poorly with the latter. This suggests that the crowding of ring substituents upon chain formation alters the ring flexibility in a nonsystematic fashion. In general, the description of carbohydrate chains afforded by 56A6(CARBO_R) suggests a significant extent of ring flexibility, i.e., small but often non-negligible equilibrium populations of inverted chairs, and challenges the "textbook" picture of conformationally locked carbohydrate rings.

European Biophysics Journal, 42 (7), ISSN:0175-7571, ISSN:1432-1017

The absolute intrinsic hydration free energy GH+,wat◦ of the proton, the surface electric potential jump χwat◦ upon entering bulk water, and the absolute redox potential VH+,wat◦ of the reference hydrogen electrode are cornerstone quantities for formulating single-ion thermodynamics on absolute scales. They can be easily calculated from each other but remain fundamentally elusive, i.e., they cannot be determined experimentally without invoking some extra-thermodynamic assumption (ETA). The Born model provides a natural framework to formulate such an assumption (Born ETA), as it automatically factors out the contribution of crossing the water surface from the hydration free energy. However, this model describes the short-range solvation inaccurately and relies on the choice of arbitrary ion-size parameters. In the present study, both shortcomings are alleviated by performing first-principle calculations of the hydration free energies of the sodium (Na+) and potassium (K+) ions. The calculations rely on thermodynamic integration based on quantum-mechanical molecular-mechanical (QM/MM) molecular dynamics (MD) simulations involving the ion and 2000 water molecules. The ion and its first hydration shell are described using a correlated ab initio method, namely resolution-of-identity second-order Moller-Plesset perturbation (RIMP2). The next hydration shells are described using the extended simple point charge water model (SPC/E). The hydration free energy is first calculated at the MM level and subsequently increased by a quantization term accounting for the transformation to a QM/MM description. It is also corrected for finite-size, approximate-electrostatics, and potential-summation errors, as well as standard-state definition. These computationally intensive simulations provide accurate first-principle estimates for GH+,wat◦, χwat◦, and VH+,wat◦, reported with statistical errors based on a confidence interval of 99%. The values obtained from the independent Na+ and K+ simulations are in excellent agreement. In particular, the difference between the two hydration free energies, which is not an elusive quantity, is 73.9 ± 5.4 kJ mol-1 (K+ minus Na+), to be compared with the experimental value of 71.7 ± 2.8 kJ mol-1. The calculated values of GH+,wat◦, χwat◦, and VH+,wat◦ (-1096.7 ± 6.1 kJ mol-1, 0.10 ± 0.10 V, and 4.32 ± 0.06 V, respectively, averaging over the two ions) are also in remarkable agreement with the values recommended by Reif and Hunenberger based on a thorough analysis of the experimental literature (-1100 ± 5 kJ mol-1, 0.13 ± 0.10 V, and 4.28 ± 0.13 V, respectively). The QM/MM MD simulations are also shown to provide an accurate description of the hydration structure, dynamics, and energetics.

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 recent article (Bieler et al., J. Chem. Theory Comput. 2014, 10, 3006), we introduced a combination of λ-dynamics and local-elevation umbrella-sampling termed λ-LEUS to calculate free-energy changes associated with alchemical processes using molecular dynamics simulations. This method was suggested to be more efficient than thermodynamic integration (TI), because the dynamical variation of the alchemical variable λ opens up pathways to circumvent barriers in the orthogonal space (defined by the N – 1 degrees of freedom that are not subjected to the sampling enhancement), a feature λ-LEUS shares with Hamiltonian replica-exchange (HR) approaches. However, the mutation considered, hydroquinone to benzene in water, was no real challenge in terms of orthogonal-space properties, which were restricted to solvent-relaxation processes. In the present article, we revisit the comparison between TI and λ-LEUS considering non-trivial mutations of the central residue X of a KXK tripeptide in water (with X = G, E, K, S, F, or Y). Side-chain interactions that may include salt bridges, hydrogen bonds or steric clashes lead to slow relaxation in the orthogonal space, mainly in the two-dimensional subspace spanned by the central φ and ψ dihedral angles of the peptide. The efficiency enhancement afforded by λ-LEUS is confirmed in this more complex test system and can be attributed explicitly to the improved sampling of the orthogonal space. The sensitivity of the results to the nontrivial choices of a mass parameter and of a thermostat coupling time for the alchemical variable is also investigated, resulting in recommended ranges of 50 to 100 u nm2 and 0.2 to 0.5 ps, respectively. © 2015 Wiley Periodicals, Inc.

The effect of methanol on the phase and phase-transition properties of a 2×8×8 glycerol-1-monopalmitate bilayer patch is investigated using a series of 239 molecular dynamics simulations on the 180 ns timescale, considering methanol concentrations cM and temperatures T in the ranges 0-12.3M and 302-338 K, respectively. The results in the form of hysteresis-corrected transition temperatures Tm are compatible with the expected features of the biphasic effect, with a reversal concentration crev of about 5.2 M. Below this concentration, the main transition is between the liquid crystal (LC) and gel (GL) phases, and Tm decreases upon increasing cM. Above this concentration, the interdigitated (ID) phase is the stable ordered phase instead, and Tm slightly increases upon increasing T up to about 10 M. The analysis of the structural and dynamical properties also reveals very different sensitivities and responses of the three phases to changes in cM. In particular, the properties of the GL phase are insensitive to cM, whereas those of the LC and ID phases are altered via an increase of the area per lipid. For the LC phase, increasing cM promotes disorder and fluidity. For the ID phase, in contrast, increasing cM up to about 10 M slightly increases the ordering and rigidity. Two side issues are also addressed, concerning: (i) the occurrence tilt-precession motions in the GL and ID phases; (ii) the influence of the pressure coupling scheme employed in the simulations, semi- or fully-anisotropic, on the simulation results.

The occurrence of long-timescale motions in glycerol-1-monopalmitate (GMP) lipid bilayers is investigated based on previously reported 600 ns molecular dynamics simulations of a 2 × 8 ×8 GMP bilayer patch in the temperature range 302–338 K, performed at three different hydration levels, or in the presence of the cosolutes methanol or trehalose at three different concentrations. The types of long-timescale motions considered are: (i) the possible phase transitions; (ii) the precession of the relative collective tilt-angle of the two leaflets in the gel phase; (iii) the trans–gauche isomerization of the dihedral angles within the lipid aliphatic tails; and (iv) the flipping of single lipids across the two leaflets. The results provide a picture of GMP bilayers involving a rich spectrum of events occurring on a wide range of timescales, from the 100-ps range isomerization of single dihedral angles, via the 100-ns range of tilt precession motions, to the multi-μs range of phase transitions and lipid-flipping events.

When using molecular dynamics (MD) simulations for the calculation of alchemical free energy changes the extended dynamics method called ? dynamics (?D) in its currently available implementations does not compare very favorably with thermodynamic integration (TI) in terms of robustness efficiency and accuracy although it is in principle easier to set up postprocess and automatize. In the present article the main shortcomings of the ?D approach are carefully analyzed and possible remedies are proposed. The resulting scheme called ? LEUS involves: (i) the use of a simple noninvertible coordinate transformation ensuring a finite sampling of the two physical end states; (ii) the application of the local elevation umbrella sampling (LEUS) memory based biasing scheme to enforce homogeneous sampling and overcome barriers along the alchemical coordinate; (iii) recommendations concerning the choice of the mass parameter and of the temperature coupling scheme for this coordinate; (iv) the use of a second order splines basis set for the memory based biasing functions. The ? LEUS scheme is described and tested considering the perturbation of hydroquinone to benzene in water. The results are compared to those of TI calculations and exhibit a superior accuracy to efficiency ratio presumably because dynamic variations in the alchemical coordinate open up pathways to circumvent orthogonal barriers these pathways being inaccessible when the coordinate is constrained. Therefore ? LEUS combines the practical advantages of ?D with the robustness of TI simultaneously affording a slightly enhanced computational efficiency. © 2014 American Chemical Society.

Phase-transition properties of glycerol-1-monopalmitate (GMP) bilayers are investigated using explicit-solvent molecular dynamics (MD) simulations, initiated from structures appropriate for the gel (GL) or liquid crystal (LC) phases, and carried out at different hydration levels and temperatures. Building up on a previous study and based on 600 ns simulations, the influence of the system size and of the force field on the equilibrium thermodynamic and dynamic parameters of the bilayers in the GL and LC phases, as well as on the temperature Tm and properties of the GL ↔ LC phase transition, are analysed. Qualitatively speaking, the results agree with the available experimental data for the area per lipid in the two phases and for the phase-transition temperatures at the three hydration levels irrespective of the selected model parameters. They also suggest that the total number of hydrogen bonds formed between a lipid headgroup and its environment is essentially constant, amounting to about four in both th...