Article pubs.acs.org/JPCB
Use of Enveloping Distribution Sampling to Evaluate Important Characteristics of Biomolecular Force Fields Wei Huang, Zhixiong Lin, and Wilfred F. van Gunsteren* Laboratory of Physical Chemistry, Swiss Federal Institute of Technology, ETH, Vladimir-Prelog-Weg 2/HCI, 8093 Zürich, Switzerland S Supporting Information *
ABSTRACT: The predictive power of biomolecular simulation critically depends on the quality of the force field or molecular model used and on the extent of conformational sampling that can be achieved. Both issues are addressed. First, it is shown that widely used force fields for simulation of proteins in aqueous solution appear to have rather different propensities to stabilize or destabilize α-, π-, and 310helical structures, which is an important feature of a biomolecular force field due to the omni-presence of such secondary structure in proteins. Second, the relative stability of secondary structure elements in proteins can only be computationally determined through so-called free-energy calculations, the accuracy of which critically depends on the extent of configurational sampling. It is shown that the method of enveloping distribution sampling is a very efficient method to extensively sample different parts of configurational space.
■
umbrella biasing MD17 requires a choice of pathway from conformation α to conformation β that should not involve high free-energy barriers, along which multiple MD simulations are to be carried out.10,14,18−21 Recently, a method called enveloping distribution sampling22−28 (EDS) has been developed that offers a viable, efficient alternative to these two standard methods to obtain ΔGβα values. In the EDS method, no pathway of intermediate states needs to be chosen, and the end-state conformations α and β can be equally sampled in a single simulation.29 It constitutes an efficient way to calculate free-enthalpy differences between different conformational states.27,30,31 Here we illustrate the efficiency of the EDS method by investigating the free-enthalpy differences ΔGβα between 310-, α-, and π-helical conformations of an alanine deca-peptide solvated in water for three different force fields, that is, the GROMOS force fields 45A3,32 53A6OXY+N,33,34 and 54A7,35 and compare the results with those previously calculated27 using the GROMOS force field 53A636 and to results for comparable peptides using other biomolecular force fields and slightly different simulation parameter settings as found in the literature.7−15 It turns out that the relative free enthalpies of the three helical conformations are quite sensitive to some forcefield parameter changes (Table 1), making these values rather important characteristics of a force field.
INTRODUCTION Most of the residues in known proteins are found in regular secondary structures, such as helices, sheets, and turns. Among these conformational motifs, the helix is the major building block of proteins.1,2 Three different types of helices have been found in proteins: the 310-, α-, and π-helices, which differ in number of residues per turn, that is, 3, 4, and 5, and thus in number of inter-residue hydrogen bonds for a poly peptide of a given length.2−4 Experimentally, it is not straightforward to determine the relative free enthalpy of these three helical conformations,5,6 yet an analysis of their occurrence in protein crystal structures suggests that the α-helix is the most stable one. This suggestion may depend on the force field used in the protein structure refinement process that converts X-ray diffraction intensities or NMR spectroscopic signal intensities into molecular structure. In computational studies using molecular dynamics (MD) simulations, different biomolecular force fields show a wide variety of propensities for 310-, α-, and π-helices due to differences in backbone force-field parameters; see refs 7−15 and references therein. This makes the freeenthalpy differences between these three helical conformations an important feature of any biomolecular force field. Computationally, the calculation of a free-enthalpy difference ΔGβα = Gβ − Gα between two different conformations or conformational states α and β of a polypeptide in aqueous solution is not trivial. Because in standard MD simulations the likelihood of sampling a high-energy conformation decreases exponentially with its energy U that is proportional to exp(−U/ kBT), where kB is Boltzmann’s constant and T is the temperature, a calculation of ΔGβα either requires very long simulations to sufficiently sample high-energy helical conformations or when using thermodynamic integration16 (TI) or © 2014 American Chemical Society
Special Issue: William C. Swope Festschrift Received: November 8, 2013 Revised: January 6, 2014 Published: January 10, 2014 6424
dx.doi.org/10.1021/jp411005x | J. Phys. Chem. B 2014, 118, 6424−6430
The Journal of Physical Chemistry B
Article
Table 1. Free-Enthalpy Differences ΔGπα=Gπ − Gα and ΔG310α=G310 − Gα of the Two Conformational Sets β (π-helix or 310helix) and α (α-helix) of the Alanine Deca-Peptide in Aqueous Solution force field
45A3
53A6
53A6OXY+N
54A7
ΔGπα [kJ/mol] ΔG310α [kJ/mol]
5.0 ± 0.5 46.7 ± 1.3
5.0 ± 0.8 47.1 ± 2.5
0.8 ± 1.1 38.5 ± 1.0
14.3 ± 1.1 36.7 ± 1.2
The GROMOS force field37 is meant for use in simulations at ambient temperatures and pressures of liquids, crystals, and solutions, in particular, in the context of biomolecular systems such as peptides and proteins, nucleic acids, carbohydrates, and lipids. The most widely used versions of the GROMOS force field are the 43A1 force field38,39 of 1996, the 45A3 force field32 of 2001, the 53A6 force field36 of 2004, and the 54A7 force field35 of 2011. The force field 43A1 contains 43 individual atom types to describe van der Waals interactions, and the van der Waals parameters of the aliphatic and aromatic carbons were obtained by calibration against thermodynamic data such as density, heat of vaporization, and free enthalpy of hydration for short alkanes up to hexane.39 The 45A3 parameter set was a reparametrization against such thermodynamic properties for a larger set of short and long alkane chains.32 In the subsequent development of the GROMOS force field, leading to version 53A6,36 the nonbonded interaction parameters of atoms in polar amino acid side chains and the peptide backbone moiety were calibrated to the mentioned type of thermodynamic data of amino-acid analogues. The 53A6OXY+N parameter set was a reparametrization for a larger set of small molecular compounds containing oxygen and nitrogen atoms.33,34 Although it is based on a sound thermodynamic fundament, the 53A6 force field turned out to insufficiently maintain α-helical structure in proteins. This induced a slight modification incorporated in the 54A735 force field.
corresponding helical hydrogen-bond populations are about 60−80% (Tables S1−S3, Supporting Information). The conformational sets α and β corresponding to the helices were defined through the atom-positional root-meansquare deviation (RMSD) of the backbone atoms (N, CA, C) of the peptide (including the two termini) from the ideal helix. Conformations belong to set α if RMSD( r ⃗ N , rα⃗ N ) ≤ RMSDαthres and RMSD( r ⃗ N , r β⃗ N ) > RMSDthres β (2)
and they belong to set β if RMSD( r ⃗ N , rα⃗ N ) > RMSDαthres and RMSD( r ⃗ N , r β⃗ N ) ≤ RMSDthres β (3)
The three ideal helices (rNξ⃗ ), with ξ = α, π, or 310, were defined through the φ and ψ backbone torsional-angle values −57.8 and −47.0° for the α-helix, −57.0 and −70.0° for the π-helix, and −49.0 and −27.0° for the 310-helix, followed by an energyminimization in vacuo using the GROMOS 53A6 force field. The resulting configurations were used as reference conformations rαN⃗ , rπN⃗ , and r3N⃗ 10. The RMSD threshold value RMSDthres was set to 0.15 nm for all three helices (ξ = α, π, ξ or 310). End-State Simulations and Enveloping Distribution Sampling Simulations. The end-state simulations for the three helices (310-, α-, and π-helices) were carried out for 11 ns using three different force fields Vphys(rN⃗ ) and the hydrogenbond restraining potential energy terms (eq 1). The potential energy of the EDS reference-state Hamiltonian VR of the N particles of the system is defined as
■
MATERIALS AND METHODS Molecular Model, Definition of End-State Hamiltonians, and Conformational Sets. The model system considered is an alanine deca-peptide capped at both termini with methyl groups, acetyl-(Ala)10-N-methyl, solvated in water. The GROMOS force fields 45A3,32 53A6OXY+N,33 and 54A735 were used for the peptide, respectively. The simple point charge (SPC) model40 was used to model water. N rest N The restraining potential energy term Vrest X (r ⃗ ; KX , r0⃗ X) used to characterize different end-state Hamiltonians is defined as an attractive half-harmonic function applied to the hydrogenbonding pairs of O and H atoms characterizing helical conformations: V Xrest( r ⃗ N ; KXrest , r0⃗ NX ) =
1 rest KX 2
=0
rest
R
R ) + V phys( r ⃗ N ) = V EDS,rest( r ⃗ N ; s , E BA
(4)
Vrest A
Vrest B
where and are the restraining potential energy functions of the two end states A and B, s is a smoothness parameter, and ERB − ERA = ERBA is an energy offset parameter difference. s and ERBA are determined such as to optimize the sampling of both end states A and B in a single simulation of VR. The EDS parameter update simulations were performed to find optimal reference state parameters for calculating the freeenthalpy difference between the π-helix and the α-helix and between the 310-helix and the α-helix for 128 × 0.15 ns. The πhelix and the α-helix served as the initial configuration, respectively. The parameters s and ERB = ERBA (ERA is standardly set to zero in two-state EDS) were updated at fixed time points: after the 1st, 3rd, 7th, 15th, and so on of 0.15 ns simulation periods. That is, the simulation time period was doubled after each update. The detailed updating scheme is described in ref 24. The resulting s and ERB parameters were used for EDS evaluation simulations of 51 ns. The free-enthalpy differences
NHB, X
∑
N
R ) = −kBTs−1 ln{e−s(V A ( r ⃗ ) − EA )/ kBT V ( r ⃗ N ; s , E BA R rest N + e−s(VB ( r ⃗ ) − EB )/ kBT } + V phys( r ⃗ N )
(diX − d0X )2 when diX > d0X
i=1
when diX ≤ d0X
(1)
where X = A, B, and C, which restrains the configuration r ⃗ of the N atoms of the peptide to the configuration rN0⃗ X of an α-, π-, or 310-helix, respectively. NHB,X is the number of hydrogen bonds (8, 7, or 9) of these three helices, and diX is the distance between the hydrogen bonding O and H atoms in the configurations in the simulation using restraining function Vrest X (Tables S1−S3, Supporting Information). The reference distance d0X was set to 0.25 nm for states A and B and 0.19 nm for state C. Krest X is the force constant set to different values to ensure that in the end-state simulations the averages of the N
6425
dx.doi.org/10.1021/jp411005x | J. Phys. Chem. B 2014, 118, 6424−6430
The Journal of Physical Chemistry B
Article
between the two conformational sets α and β, that is, ΔGβα = Gβ − Gα, were calculated through24 ⎧ +V EDS,rest / kBT >R ,set β Nβ(VR ) ⎫ ⎪ R ,set α Nα(VR ) ⎭ ⎩
