An effective all-atom potential for proteins.

We describe and test an implicit solvent all-atom potential for simulations of protein folding and aggregation. The potential is developed through studies of structural and thermodynamic properties of 17 peptides with diverse secondary structure. Results obtained using the final form of the potential are presented for all these peptides. The same model, with unchanged parameters, is furthermore applied to a heterodimeric coiled-coil system, a mixed alpha/beta protein and a three-helix-bundle protein, with very good results. The computational efficiency of the potential makes it possible to investigate the free-energy landscape of these 49-67-residue systems with high statistical accuracy, using only modest computational resources by today's standards.PACS Codes: 87.14.E-, 87.15.A-, 87.15.Cc.

1 Introduction

A molecular understanding of living systems requires modeling of the dynamics and interactions of n class="Chemical">proteins. The relevant dynamics of a protein may amount to small fluctuations about its native structure, or reorientations of its ordered parts relative to each other. In either case, a tiny fraction of the conformational space is explored. For flexible proteins, perhaps with large intrinsically disordered parts [1,2], the situation is different. When studying such proteins or conformational conversion processes like folding or amyloid aggregation, the competition between different minima on the free-energy landscape inevitably comes into focus. Studying these systems by computer simulation is a challenge, because proper sampling of all relevant free-energy minima must be ensured. This goal is very hard to achieve if explicit solvent molecules are included in the simulations. The use of coarse-grained models can alleviate this problem, but makes important geometric properties like secondary structure formation more difficult to describe.

Here we present an implicit solvent all-atom protein model especially aimed at n class="Chemical">problems requiring exploration of the global free-energy landscape. It is based on a computationally convenient effective potential, with parameters determined through full-scale thermodynamic simulations of a set of experimentally well characterized peptides. Central to the approach is the use of a single set of model parameters, independent of the protein studied. This constraint is a simple but efficient way to avoid unphysical biases, for example, toward either α-helical or β-sheet structure [3,4]. Imposing this constraint is also a way to enable systematic refinement of the potential.

An earlier version [5,6] of tn class="Chemical">his potential has proven useful, for example, for studies of aggregation [7-9] and mechanical unfolding [10,11]. Also, using a slightly modified form of the potential [12], the folding mechanisms of a 49-residue protein, Top7-CFr, were investigated [13,14]. Here we revise this potential, through studies of an enlarged set of 17 peptides (see Table 1 and Fig. 1). We show that the model, in its final form, folds these different sequences to structures similar to their experimental structures, using a single set of potential parameters. The description of each peptide is kept brief, to be able to discuss all systems and thereby address the issue of transferability in a direct manner. The main purpose of this study is model development rather than detailed characterization of individual systems.

Table 1

Amino acid sequences

System	PDB code	Sequence
Trp-cage	1L2Y	NLYIQ WLKDG GPSSG RPPPS
E6apn1	1RIJ	Ac-ALQEL LGQWL KDGGP SSGRP PPS-NH2
C		Ac-KETAA AKFER AHA-NH2
EK		Ac-YAEAA KAAEA AKAF-NH2
Fs		Suc-AAAAA AAARA AAARA AAARA A-NH2
GCN4tp	2OVN	NYHLE NEVAR LKKLV GE
HPLC-6	1WFA	DTASD AAAAA ALTAA NAKAA AELTA ANAAA AAAAT AR-NH2
Chignolin	1UAO	GYDPE TGTWG
MBH12	1J4M	RGKWT YNGIT YEGR
GB1p		GEWTY DDATK TFTVT E
GB1m2		GEWTY NPATG KFTVT E
GB1m3		KKWTY NPATG KFTVQ E
trpzip1	1LE0	SWTWE GNKWT WK-NH2
trpzip2	1LE1	SWTWE NGKWT WK-NH2
betanova		RGWSV QNGKY TNNGK TTEGR
LLM		RGWSL QNGKY TLNGK TMEGR
beta3s		TWIQN GSTKW YQNGS TKIYT
AB zipper	1U2U	Ac-EVAQL EKEVA QLEAE NYQLE QEVAQ LEHEG-NH2
		Ac-EVQAL KKRVQ ALKAR NYALK QKVQA LRHKG-NH2
Top7-CFR	2GJH	ERVRI SITAR TKKEA EKFAA ILIKV FAELG YNDIN VTWDG DTVTV EGQL
GS-α3 W	1LQ7	GSRVK ALEEK VKALE EKVKA LGGGG RIEEL KKKWE ELKKK IEELG GGGEV KKVEE EVKKL EEEIK KL

Suc stands for succinylic acid.

Figure 1

Schematic illustration of native geometries studied. (a) the Trp-cage, (b) an α-helix, (c) a β-hairpin, (d) a three-stranded β-sheet, (e) an α-helix dimer (1U2U), (f) a three-helix bundle (1LQ7), and (g) a mixed α/β protein (2GJH).

Schematic illustration of native geometries studied. (a) the Trp-cage, (b) an α-helix, (c) a β-hairpin, (d) a n class="Chemical">three-stranded β-sheet, (e) an α-helix dimer (1U2U), (f) a three-helix bundle (1LQ7), and (g) a mixed α/β protein (2GJH).

Whether or not this potential, calibrated un class="Chemical">sing data on peptides with typically ~20 residues, will be useful for larger systems is not obvious. Therefore, we also apply our potential, with unchanged parameters, to three larger systems with different geometries. These systems are the mixed α/β protein Top7-CFr, a three-helix-bundle protein with 67 residues, and a heterodimeric leucine zipper composed of two 30-residue chains.

Protein folding n class="Chemical">simulations are by necessity based on potentials whose terms are interdependent and dependent on the choice of geometric representation. Therefore, we choose to calibrate our potential directly against folding properties of whole chains. To make this feasible, we deliberately omit many details included in force fields like Amber, CHARMM and OPLS (for a review, see [15]). With this approach, we might lose details of a given free-energy minimum, but, by construction, we optimize the balance between competing minima.

Two potentials somewhat sin class="Chemical">milar in form to ours are the μ-potential of the Shakhnovich group [16] and the PFF potential of the Wenzel group [17]. These groups also consider properties of entire chains for calibration, but use folded PDB structures or sets of decoys rather than full-scale thermodynamic simulations. Our admittedly time-consuming procedure implies that our model is trained on completely general structures, which might be an advantage when studying the dynamics of folding. Another potential with similarities to ours is that developed by the Dokholyan group for discrete molecular dynamics simulations [18].

2 Methods

Our model belongs to the class of implicit solvent all-atom models with torsional degrees of freedom. All geometrical parameters, like bond lengths and bond angles, are as dn class="CellLine">escribed earlier [5].

The interaction potential is composed of four major terms:

The first term, Eloc, contains local interactions between class="Disease">n atoms separated by only a few covalent bonds. The other three terms are non-local in character: Eev represents excluded-volume effects, Ehb is a hydrogen-bond potential, and Esc contains residue-specific interactions between pairs of sidechains. Next we describe the precise form of these four terms. Energy parameters are given in a unit called eu. The factor for conversion from eu to kcal/mol will be determined in the next section, by calibration against the experimental melting temperature for one of the peptides studied, the Trp-cage.

2.1 Local potential

The local potential can be divided into two backbone terms, and , and one sidechain term, . In dn class="CellLine">escribing the potential, the concept of a peptide unit is useful. A peptide unit consists of the backbone C'O group of one residue and the backbone NH group of the next residue.

• The potential represents interactions between partial charges of neighboring peptide units along the chain. It is given by

where the outer sum runs over all pairs of nearest-neighbor peptide units and each of the two inner sums runs over atoms in one peptide unit (if the N n class="Chemical">side of the peptide unit is proline the sum runs over only C' and O). The partial charge qis taken as ± 0.42 for C' and O atoms and ± 0.20 for H and N atoms. The parameter is set to 6 eu, corresponding to a dielectric constant of ϵr ≈ 41. Two peptide units that are not nearest neighbors along the chain interact through hydrogen bonding (see below) rather than through the potential .

• The term provides an additional OO and HH repulsion for neighboring peptide units, unless the residue flanked by the two peptide units is a glycine. This repulsion is added to make doubling of hydrogen bonds less likely. Glycine has markedly different backbone energetics compared to other residues. The lack of Catom makes glycine more flexible. However, the observed distribution of Ramachandran ϕ, ψ angles for glycine in PDB structures [19] is not as broad as simple steric considerations would suggest. provides an energy penalty for glycine ψ values around ± 120°, which are sterically allowed but relatively rare in PDB structures.

The full expression for is

where = 1.2 eu, = -0.15 eu, I is a residue index, and

The function f(u) is positive if the HHdistance, d(H, H), is smaller than both of the Hn class="Chemical">Nand NHdistances, and zero otherwise. This term thus provides an energy penalty when Hand Hare exposed to each other (it is omitted if residue I or I + 1 is a proline). Similarly, f(v) is positive when Oand Oare exposed to each other.

• is an explicit torsion angle potential for n class="Chemical">sidechain angles, χ. Many sidechain angles display distributions resembling what one would expect based on simple steric considerations. The use of the torsion potential is particularly relevant for χ2 in asparagine and aspartic acid and χ3 in glutamine and glutamic acid. The torsion potential is defined as

where and nare constants. Each sidechain angle χbelongs to one of four classes associated with different n class="Chemical">values of and n(see Table 2).

Table 2

Classification of sidechain angles, χ

Residue	χ1	χ2	χ3	χ4
Ser, Cys, Thr, Val	I
Ile, Leu	I	I
Asp, Asn	I	IV
His, Phe, Tyr, Trp	I	III
Met	I	I	II
Glu, Gln	I	I	IV
Lys	I	I	I	I
Arg	I	I	I	III

The parameters of the torsion angle potential are (, n) = (0.6 eu, 3) for class I, (, n) = (0.3 eu, 3) for class II, (, n) = (0.4 eu, 2) for class III, and (, n) = (-0.4 eu, 2) for class IV.

Classification of n class="Chemical">sidechain angles, χ

2.2 Excluded volume

Excluded-volume effects are modeled using the potential

where the summation is over all pairs of atoms with a non-constant separation, κev = 0.10 eu, and σ= 1.77, 1.75, 1.53, 1.42 and 1.00 Å for S, C, N, O and n class="Disease">H atoms, respectively. The parameter λis unity for pairs connected by three covalent bonds and λ= 0.75 for all other pairs. To speed up the calculations, Eev is evaluated using a cutoff of 4.3 λÅ.

2.3 Hydrogen bonding

Our potential contains an explicit hydrogen-bond term, Ehb. All n class="Chemical">hydrogen bonds in the model are between NH and CO groups. They connect either two backbone groups or a charged sidechain (aspartic acid, glutamic acid, lysine, arginine) with a backbone group. Two neighboring peptide units, which interact through the local potential (see above), are not allowed to hydrogen bond with each other.

The form of the hydrogen-bond potential is

where = 3.0 eu and = 2.3 eu set the strengths of backbone-backbone and sidechain-backbone bonds, respectively, ris the HO distance, αis the n class="Disease">NHO angle, and βis the HOC angle. The functions u(r) and v(α, β) are given by

where σhb = 2.0 Å. A 4.5 Å cutoff is used for u(r).

2.4 Sidechain potential

Our sidechain potential is composed of two terms, n class="CellLine">Esc = Ehp + Ech. The Ech term represents interactions among sidechain charges. The first and more important term, Ehp, is meant to capture the effects of all other relevant interactions, especially effective hydrophobic attraction. For convenience, Ehp and Ech have a similar form,

Here the sums run over residue pairs IJ, and are contact measures that take n class="Chemical">values between 0 and 1, and and are energy parameters.

It is assumed that ten of the twenty natural amino acids contribute to Ehp, see Table 3. Included among these ten are n class="Chemical">lysine and arginine, which are charged but have large hydrophobic parts. To reduce the number of parameters, the hydrophobic contact energies are taken to be additive, = m+ m. It is known that the statistically derived Miyazawa-Jernigan contact matrix [20] can be approximately decomposed this way [21]. The mparameters can be found in Table 3. is set to 0 if residues I and J are nearest neighbors along the chain, and is reduced by a factor 2 for next-nearest neighbors.

Table 3

The parameter mof the hydrophobicity potential Ehp

Residue	mI (eu)
Arg	0.3
Met, Lys	0.4
Val	0.6
Ile, Leu, Pro	0.8
Tyr	1.1
Phe, Trp	1.6

The residues taken as chn class="Chemical">arged are aspartic acid, glutamic acid, lysine and arginine. The charge-charge contact energy is – = 1.5sseu, where sand sare the signs of the charges (± 1).

The contact measure is calculated using a predetern class="Chemical">mined set of atoms for each amino acid, denoted by (see Table 4). Let nbe the number of atoms in and let

Table 4

Atoms used in the calculation of the contact measure

Residue	Set of atoms (AI)
Pro	Cβ, Cγ, Cδ
Tyr	Cγ, Cδ1, Cδ2, Cϵ1, Cϵ2, Cζ
Val	Cβ, Cγ1, Cγ2
Ile	Cβ, Cγ1, Cγ2, Cδ
Leu	Cβ, Cγ, Cδ1, Cδ2
Met	Cβ, Cγ, Sδ, Cϵ
Phe	Cγ, Cδ1, Cδ2, Cϵ1, Cϵ2, Cζ
Trp	Cγ, Cδ1, Cδ2, Cϵ3, Cζ3, Cη2
Arg	Cβ, Cγ
Lys	Cβ, Cγ, Cδ

where g(x) is unity for x < (3.7 Å)2, vanishes for x > (4.5 Å)2, and varies linearly for intermediate x. The contact measure can then be written as

where γis either 1 or 0.75. For γ= 1, is, roughly speaking, the fraction of atoms in and that are in contact with some atom from the other of the two sets. A reduction to γ= 0.75 makes it easier to achieve a full contact ( = 1). The n class="Chemical">value γ= 0.75 is used for interactions within the group proline, phenylalanine, tyrosine and tryptophan, to make face-to-face stacking of these sidechains less likely. It is also used within the group isoleucine, leucine and valine, because a full contact is otherwise hard to achieve for these pairs. In all other cases, γis unity.

The definition of is sin class="Chemical">milar. The γparameter is unity for charge-charge interactions, and the sets of atoms used, , can be found in Table 5.

Table 5

Atoms used in the calculation of the contact measure

Residue	Set of atoms (AI)
Arg	Nϵ, Cζ, Nη1, Nη2
Lys	1Hζ, 2Hζ, 3Hζ
Asp	Oδ1, Oδ2
Glu	Oϵ1, Oϵ2

2.5 Chain ends

Some of the sequences we study have extra groups attached at one or both ends of the chain. The groups occurring are N-tern class="Chemical">minal acetyl and succinylic acid, and C-terminal NH2. When such a unit is present, the model assumes polar NH and CO groups beyond the last Catom to hydrogen bond like backbone NH/CO groups but with the strength reduced by a factor 2 (multiplicatively). The charged group of succinylic acid interacts like a charged sidechain.

In the absence of end groups, the model assumes the N and C tern class="Chemical">mini to be positively and negatively charged, respectively, and to interact like charged sidechains.

2.6 Monte Carlo details

We investigate the folding thermodynamics of tn class="Chemical">his model by Monte Carlo (MC) methods. The simulations are done using either simulated tempering (ST) [22,23] or parallel tempering/replica exchange (PT) [24,25], both with temperature as a dynamical variable. For small systems we use ST, with seven geometrically distributed temperatures in the range 279 K–367 K. For each system, ten independent ST runs are performed. For our largest systems we use PT with a set of sixteen temperatures, spanning the same interval. Using fourfold multiplexing [26], one run comprising 64 parallel trajectories is performed for each system. The PT temperature distribution is determined by an optimization procedure [26]. The length of our different simulations can be found in Table 6.

Table 6

Algorithm used and total number of elementary MC steps for all systems studied

System	Method	MC steps
Trp-cage, E6apn1	ST	10 × 1.0 × 109
C, EK, Fs, GCN4tp	ST	10 × 1.0 × 109
HPLC-6	ST	10 × 3.0 × 109
Chignolin	ST	10 × 0.5 × 109
MBH12	ST	10 × 1.0 × 109
GB1p	ST	10 × 2.0 × 109
GB1m2, GB1m3	ST	10 × 1.0 × 109
Trpzip1, trpzip2	ST	10 × 1.0 × 109
betanova, LLM	ST	10 × 1.0 × 109
beta3s	ST	10 × 2.0 × 109
AB zipper	PT	64 × 3.0 × 109
Top7-CFR	PT	64 × 2.4 × 109
GS-α3 W	PT	64 × 3.5 × 109

Three different conformational updates are used in the n class="Chemical">simulations: single variable updates of sidechain and backbone angles, respectively, and Biased Gaussian Steps (BGS) [27]. The BGS move is semi-local and updates up to eight consecutive backbone degrees of freedom in a manner that keeps the ends of the segment approximately fixed. The ratio of sidechain to backbone updates is the same at all temperatures, whereas the relative frequency of the two backbone updates depends on the temperature. At high temperatures the single variable update is the only backbone update used, and at low temperatures only BGS is used. At intermediate temperatures both updates are used.

The AB zipper, a two-chain system, is studied using a periodic box of n class="Chemical">size (158 Å)3. In addition to the conformational updates described above, the simulations of this system used rigid body translations and rotations of individual chains.

Our simulations are performed un class="Chemical">sing the open source C++-package PROFASI [28]. Future public releases of PROFASI will include an implementation of the force field described here. While this force field has been implemented in PROFASI in an optimized manner, this optimization does not involve a parallel evaluation of the potential on many processors. Therefore, in our simulations the number of processors used is the same as the number of MC trajectories generated. For a typical small peptide, a trajectory of the length as given in Table 6 takes ~18 hours to generate on an AMD Opteron processor with ~2.0 GHz clock rate. For the largest system studied, GS-α3 W, the simulations, with a proportionately larger number of MC updates, take ~10 days to complete.

2.7 Analysis

In our simulations, we monitor a variety of different n class="Chemical">properties. Three important observables are as follows.

1. α-helix content, h. A residue is defined as helical if its Ramachandran angle pair is in the region -90° <ϕ < -30°, -77° <ψ < -17°. Following [29], a stretch of n > 2 helical ren class="Chemical">sidues is said to form a helical segment of length n - 2. For an end residue that is not followed by an extra end group, the (ϕ, ψ) pair is poorly defined. Thus, for a chain with N residues, the maximum length of a helical segment is N - 4, N - 3 or N - 2, depending on whether there are zero, one or two end groups. The α-helix content h is defined as the total length of all helical segments divided by this maximum length.

2. Root-mean-square deviation from a folded reference structure, bRMSD/RMSD/pRMSD. bRMSD is calculated over backbone atoms, whereas RMSD is calculated over all heavy atoms. All residues except the two end ren class="Chemical">sidues are included in the calculation, unless otherwise stated. For the case of the dimeric AB zipper, the periodic box used for the simulations has to be taken into account. The two chains in the simulation might superficially appear to be far away when they are in fact close, because of periodicity. For this case we evaluate backbone RMSD over atoms taken from both chains in the dimer, and minimize this value with respect to periodic translations. We denote this as pRMSD.

3. Nativeness measure based on n class="Chemical">hydrogen bonds, qhb. This observable has the value 1 if at most two native backbone-backbone hydrogen bonds are missing, and is 0 otherwise. A hydrogen bond is considered formed if its energy is less than -1.03 eu.

In many cases, it turns out that the temperature dependence of our results can be approximately dn class="CellLine">escribed in terms of the simple two-state model

where X(T) is the quantity studied, X1 and X2 are the values of X in the two states, and K(T) is the effective equilibrium constant (R is the gas constant). In tn class="Chemical">his first-order form, K(T) contains two parameters: the melting temperature Tm and the energy difference ΔE. The parameters Tm, ΔE, X1 and X2 are determined by fitting to data.

Thermal averages and their statistical errors are calculated by using the jackknife method [30], after discarding the first 20% of each MC trajectory for thermalization.

Figures of 3D structures were prepared using PyMOL [31].

3 Results

We study a total of 20 peptide/protein systems, listed in Table 1 (an class="Chemical">mino acid sequences can be found in this table). Among these, there are 17 smaller systems with 10–37 residues and 3 larger ones with ≥ 49 residues. Many of the smaller systems have been simulated by other groups, in some cases with explicit water (for a review, see [32]). Two of the three larger systems, as far as we know, have not been studied using other force fields. A study of the 67-residue three-helix-bundle protein GS-α3 W using the ECEPP/3 force field was recently reported [33]. The simulations presented here use the same geometric representation and find about a hundred times the number of independent folding events, while consuming much smaller computing resources.

3.1 Trp-cage and E6apn1

The Trp-cage is a den class="Chemical">signed 20-residue miniprotein with a compact helical structure [34]. Its NMR-derived native structure (see Fig. 1) contains an α-helix and a single turn of 310-helix [34]. The E6apn1 peptide was designed using the Trp-cage motif as a scaffold, to inhibit the E6 protein of papillomavirus [35]. E6apn1 is three residues larger than the Trp-cage but has a similar structure, except that the α-helix is slightly longer [35].

As indicated earlier, we use melting data for the Trp-cage to set the energy scale of the model. For tn class="Chemical">his peptide, several experiments found a similar melting temperature, Tm ~315 K [34,36,37]. In our model, the heat capacity of the Trp-cage displays a maximum at RT = 0.4722 ± 0.0008 eu. Our energy unit eu is converted to kcal/mol by setting this temperature equal to the experimental melting temperature (315 K). Having done that, there is no free parameter left in the model. Other systems are thus studied without tuning any model parameter. For E6apn1, the experimental melting temperature is Tm ~305 K [35].

Fig. 2a shows the helix content h against temperature for the Trp-cage and E6apn1, as obtained from our n class="Chemical">simulations. In both cases, the T dependence is well described by the simple two-state model of Eq. 15. The fitted melting temperatures are Tm = 309.6 ± 0.7 K and Tm = 304.0 ± 0.5 K for the Trp-cage and E6apn1, respectively. This Tm value for the Trp-cage is slightly lower than that we obtain from heat capacity data, 315 K. A fit to our data for the hydrophobicity energy Ehp (not shown) gives instead a slightly larger Tm, 321.1 ± 0.8 K. This probe dependence of Tm implies an uncertainty in the determination of the energy scale. By using the Trp-cage, this uncertainty is kept small (~2%). For many other peptides, the spread in Tm is much larger (see below).

Figure 2

The Trp-cage and E6apn1. (a) Helix content h against temperature. The lines are two-state fits (Tm = 309.6 ± 0.7 K and ΔE = 11.3 ± 0.3 kcal/mol for the Trp-cage; Tm = 304.0 ± 0.5 K and ΔE = 14.2 ± 0.3 kcal/mol for E6apn1). (b) Free energy F calculated as a function of bRMSD at two different temperatures, 279 K (solid lines) and 306 K (dashed lines). The double lines indicate the statistical errors.

The Trp-cage and E6apn1. (a) Helix content h against temperature. The lines are two-state n class="Disease">fits (Tm = 309.6 ± 0.7 K and ΔE = 11.3 ± 0.3 kcal/mol for the Trp-cage; Tm = 304.0 ± 0.5 K and ΔE = 14.2 ± 0.3 kcal/mol for E6apn1). (b) Free energy F calculated as a function of bRMSD at two different temperatures, 279 K (solid lines) and 306 K (dashed lines). The double lines indicate the statistical errors.

Fig. 2b shows the free energy calculated as a function of bRMSD for the Trp-cage and E6apn1 at two different temperatures. The first temperature, 279 K, is well below Tm. Here native-like conformations don class="Chemical">minate and the global free-energy minima are at 2.4 Å and 2.0 Å for the Trp-cage and E6apn1, respectively. At the second temperature, 306 K, the minima are shifted to higher bRMSD. Note that these free-energy profiles, taken near Tm, show no sign of a double-well structure. Hence, these peptides do not show a genuine two-state behavior in our simulations, even though the melting curves (Fig. 2a) are well described by a two-state model, as are many experimentally observed melting curves.

3.2 The α-helices C, EK, Fs, GCN4tp and HPLC-6

Our next five sequences form α-helices. Among these, there are large differences in helix stability, according to n class="Chemical">CD studies. The least stable are the C [38] and EK [39] peptides, which are only partially stable at T ~273 K. The original C peptide is a 13-residue fragment of ribonuclease A, but the C peptide here is an analogue with two alanine substitutions and a slightly increased helix stability [40]. The EK peptide is a designed alanine-based peptide with 14 residues.

Our third α-helix peptide is the 21-residue n class="Chemical">Fs [41], which is also alanine-based. Fs is more stable than C and EK [41,42], with estimated Tm values of 308 K [42] and 303 K [43] from CD studies and 334 K from an IR study [44]. Even more stable is HPLC-6, a winter flounder antifreeze peptide with 37 residues. CD data suggest that the helix content of HPLC-6 remains non-negligible, ~0.10, at temperatures as high as ~343 K [45]. Our fifth helix-forming sequence, which we call GCN4tp, has 17 residues and is taken from a study of GCN4 coiled-coil formation [46]. Its melting behavior has not been studied, as far as we know, but its structure was characterized by NMR [46].

These five peptides are indeed α-helical in our model. At 279 K, the calculated helix content h is 0.28 for the C peptide, 0.47 for the EK peptide, and > 0.60 for the other n class="Chemical">three peptides. Fig. 3 shows the temperature dependence of h. By fitting Eq. 15 to the data for the three stable sequences, we find melting temperatures of 298.9 ± 0.1 K, 309.2 ± 0.3 K and 323.3 ± 1.2 K for GNC4tp, Fs and HPLC-6, respectively.

Figure 3

The C, EK, F. Helix content h against temperature. The lines are two-state fits (Tm = 276.3 ± 2.4 K and ΔE = 11.7 ± 0.4 kcal/mol for C; Tm = 293.9 ± 0.4 K and ΔE = 12.6 ± 0.2 kcal/mol for EK; Tm = 309.2 ± 0.3 K and ΔE = 18.7 ± 0.4 kcal/mol for Fs; Tm = 298.9 ± 0.1 K and ΔE = 14.1 ± 0.1 kcal/mol for GCN4tp; Tm = 323.3 ± 1.2 K and ΔE = 23.6 ± 2.2 kcal/mol for HPLC-6).

The C, EK, F. Helix content h against temperature. The lines are two-state fits (Tm = 276.3 ± 2.4 K and ΔE = 11.7 ± 0.4 kcal/mol for C; Tm = 293.9 ± 0.4 K and ΔE = 12.6 ± 0.2 kcal/mol for EK; Tm = 309.2 ± 0.3 K and ΔE = 18.7 ± 0.4 kcal/mol for n class="Chemical">Fs; Tm = 298.9 ± 0.1 K and ΔE = 14.1 ± 0.1 kcal/mol for GCN4tp; Tm = 323.3 ± 1.2 K and ΔE = 23.6 ± 2.2 kcal/mol for HPLC-6).

For the four peptides whose melting behavior has been studied experimentally, these results are in good agreement with experimental data. In particular, we find that HPLC-6 indeed is more stable than n class="Chemical">Fs in the model, which in turn is more stable than both C and EK. The model thus captures the stability order among these peptides.

3.3 The β-hairpins chignolin and MBH12

We now turn to β-sheet peptides and begin with the β-hairpins chignolin [47] and MBH12 [48] with 10 and 14 ren class="Chemical">sidues, respectively. Both are designed and have been characterized by NMR. For chignolin, Tm values in the range 311–315 K were reported [47], based on CD and NMR. We are not aware of any melting data for MBH12.

Fig. 4 shows the temperature dependence of the hydrophobicity energy Ehp and the nativeness parameter qhb for these n class="Chemical">peptides. By fitting to Ehp data, we obtain Tm = 311.0 ± 0.5 K and Tm = 315.4 ± 1.3 K for chignolin and MBH12, respectively. Using qhb data instead, we find Tm = 305.4 ± 0.5 K for chignolin and Tm = 309.2 ± 0.7 K for MBH12. These Tm values show a significant but relatively weak probe dependence. The values for chignolin can be compared with experimental data, and the agreement is good.

Figure 4

Chignolin and MBH12. (a) Hydrophobicity energy Ehp against temperature. The lines are two-state fits (Tm = 311.0 ± 0.5 K and ΔE = 9.6 ± 0.2 kcal/mol for chignolin; Tm = 315.4 ± 1.3 K and ΔE = 9.9 ± 0.9 kcal/mol for MBH12). (b) Nativeness qhb against temperature. The lines are two-state fits (Tm = 305.4 ± 0.5 K and ΔE = 10.4 ± 0.1 kcal/mol for chignolin; Tm = 309.2 ± 0.7 K and ΔE = 13.5 ± 0.2 kcal/mol for MBH12).

Chignolin and MBH12. (a) Hydrophobicity energy Ehp against temperature. The lines are two-state fits (Tm = 311.0 ± 0.5 K and ΔE = 9.6 ± 0.2 kcal/mol for chignolin; Tm = 315.4 ± 1.3 K and ΔE = 9.9 ± 0.9 kcal/mol for MBH12). (b) n class="Chemical">Nativeness qhb against temperature. The lines are two-state fits (Tm = 305.4 ± 0.5 K and ΔE = 10.4 ± 0.1 kcal/mol for chignolin; Tm = 309.2 ± 0.7 K and ΔE = 13.5 ± 0.2 kcal/mol for MBH12).

Because these peptides have only four native n class="Chemical">hydrogen bonds each, one may question our definition of qhb (see Methods), which takes a conformation as native-like (qhb = 1) even if two hydrogen bonds are missing. Therefore, we repeated the analysis using the stricter criterion that native-like conformations (qhb = 1) may lack at most one hydrogen bond. The resulting decrease in native population, as measured by the average qhb, was ~0.1 or smaller at all temperatures. Even with this stricter definition, we find native populations well above 0.5 at low temperatures for both peptides.

3.4 The β-hairpins GB1p, GB1m2 and GB1m3

GB1p is the second β-hairpin of the B1 domain of n class="Chemical">protein G (residues 41–56). Its folded population has been estimated by CD/NMR to be 0.42 at 278 K [49] and ~0.30 at 298 K [50], whereas a Trp fluorescence study found a Tm of 297 K [51], corresponding to a somewhat higher folded population. GB1m2 and GB1m3 are two mutants of GB1p with significantly enhanced stability [50]. At 298 K, the folded population was found to be 0.74 ± 0.05 for GB1m2 and 0.86 ± 0.03 for GB1m3, based on CD and NMR measurements [50]. It was further estimated that Tm = 320 ± 2 K for GB1m2 and Tm = 333 ± 2 K for GB1m3 [50].

All these three n class="Chemical">peptides are believed to adopt a structure similar to that GB1p has as part of the protein G B1 domain (PDB code 1GB1). This part of the full protein contains seven backbone-backbone hydrogen bonds. These hydrogen bonds are the ones we consider when evaluating qhb for these peptides.

Fig. 5 shows the observables Ehp and n class="Chemical">qhb against temperature for these peptides. Fits to the data give Ehp-based Tm values of 301.7 ± 3.3 K, 324.4 ± 1.1 K and 331.4 ± 0.7 K for GB1p, GB1m2 and GB1m3, respectively, and qhb-based Tm values of 307.5 ± 0.5 K and 313.9 ± 1.4 K for GB1m2 and GB1m3, respectively. The qhb data do not permit a reliable fit for the less stable GB1p. At 298 K, we find qhb-based folded populations of 0.20, 0.64 and 0.74 for GB1p, GB1m2 and GB1m3, respectively, which can be compared with the above-mentioned experimental results (0.30, 0.74 and 0.86).

Figure 5

GB1p, GB1m2 and GB1m3. (a) Hydrophobicity energy Ehp against temperature. The lines are two-state fits (Tm = 301.7 ± 3.3 K and ΔE = 11.3 ± 1.1 kcal/mol for GB1p; Tm = 324.4 ± 1.4 K and ΔE = 13.2 ± 1.0 kcal/mol for GB1m2; Tm = 331.4 ± 0.7 K and ΔE = 14.8 ± 0.5 kcal/mol for GB1m3). (b) Nativeness qhb against temperature. The lines are two-state fits (Tm = 307.5 ± 0.5 K and ΔE = 20.7 ± 0.5 kcal/mol for GB1m2; Tm = 313.9 ± 1.4 K and ΔE = 21.4 ± 1.1 kcal/mol for GB1m3).

These results show that, in the model, the apparent folded populations of these peptides depend quite strongly on the obn class="Chemical">servable studied. Our Ehp-based results agree quite well with experimental data, especially for GB1m2 and GB1m3, whereas our qhb results consistently give lower folded populations for all peptides. The stability order is the same independent of which of the two observables we study, namely GB1p < GB1m2 < GB1m3, which is the experimentally observed order.

The stability difference between GB1m2 and n class="Gene">GB1m3 is mainly due to charge-charge interactions. In our previous model [6], these interactions were ignored, and both peptides had similar stabilities. The present model splits this degeneracy. Moreover, the magnitude of the splitting, which sensitively depends on the strength of the charge-charge interactions, is consistent with experimental data.

3.5 The β-hairpins trpzip1 and trpzip2

The 12-residue n class="Chemical">trpzip1 and trpzip2 are designed β-hairpins, each containing two tryptophans per β-strand [52]. The only difference between the two sequences is a transposition of an aspargine and a glycine in the hairpin turn. CD measurements suggest that trpzip1 and trpzip2 are remarkably stable for their size, with Tm values of 323 K and 345 K, respectively [52]. A complementary trpzip2 study, using both experimental and computational methods, found Tm values to be strongly probe-dependent [53].

Fig. 6 shows our melting curves for these peptides, based on the obn class="Chemical">servables Ehp and qhb. The Ehp-based Tm values are 319.7 ± 0.2 K and 327.1 ± 0.8 K for trpzip1 and trpzip2, respectively. Using qhb data instead, we find Tm = 303.2 ± 1.1 K for trpzip1 and Tm = 305.0 ± 1.1 K for trpzip2.

Figure 6

Trpzip1 and trpzip2. (a) Hydrophobicity energy Ehp against temperature. The lines are two-state fits (Tm = 319.7 ± 0.2 K and ΔE = 7.9 ± 0.1 kcal/mol for trpzip1; Tm = 327.1 ± 0.8 K and ΔE = 8.3 ± 0.4 kcal/mol for trpzip2). (b) Nativeness qhb against temperature. The lines are two-state fits (Tm = 303.2 ± 1.8 K and ΔE = 14.1 ± 0.5 kcal/mol for trpzip1; Tm = 305.0 ± 1.1 K and ΔE = 12.6 ± 0.3 kcal/mol for trpzip2).

Trpzip1 and n class="Chemical">trpzip2. (a) Hydrophobicity energy Ehp against temperature. The lines are two-state fits (Tm = 319.7 ± 0.2 K and ΔE = 7.9 ± 0.1 kcal/mol for trpzip1; Tm = 327.1 ± 0.8 K and ΔE = 8.3 ± 0.4 kcal/mol for trpzip2). (b) Nativeness qhb against temperature. The lines are two-state fits (Tm = 303.2 ± 1.8 K and ΔE = 14.1 ± 0.5 kcal/mol for trpzip1; Tm = 305.0 ± 1.1 K and ΔE = 12.6 ± 0.3 kcal/mol for trpzip2).

Like for the other β-hairpins discussed earlier, our qhb-based folded populations are low compared to estimates based on n class="Chemical">CD data, whereas those based on Ehp are much closer to experimental data. For trpzip2, the agreement is not perfect but acceptable, given that Tm has been found to be strongly probe-dependent for this peptide [53].

3.6 Three-stranded β-sheets: betanova, LLM and beta3s

Betanova [54], the betanova triple mutant LLM [55] and beta3s [56] are den class="Chemical">signed 20-residue peptides forming three-stranded β-sheets. All the three peptides are marginally stable. NMR studies suggest that the folded population at 283 K is 0.09 for betanova [55], 0.36 for LLM [55], and 0.13–0.31 for beta3s [56].

Fig. 7 shows our Ehp and qhb data for these n class="Chemical">peptides. From the qhb data, Tm values cannot be extracted, because the stability of the peptides is too low. At 283 K, the qhb-based folded populations are 0.08, 0.47, 0.28 for betanova, LLM and beta3s, respectively, in good agreement with the experimental results. Fits to Ehp data can be performed. The obtained Tm values are 318.8 ± 2.5 K, 305.6 ± 1.7 K and 295.7 ± 3.1 K for betanova, LLM and beta3s, respectively.

Figure 7

Betanova, LLM and beta3s. (a) Hydrophobicity energy Ehp against temperature. The lines are two-state fits (Tm = 318.8 ± 2.5 K and ΔE = 13.3 ± 2.1 kcal/mol for betanova; Tm = 305.6 ± 1.7 K and ΔE = 13.4 ± 1.0 kcal/mol for LLM; Tm = 295.7 ± 3.1 K and ΔE = 9.7 ± 0.5 kcal/mol for beta3s). (b) Nativeness qhb against temperature. Two-state fits were not possible.

Betanova, LLM and beta3s. (a) Hydrophobicity energy Ehp against temperature. The lines are two-state n class="Disease">fits (Tm = 318.8 ± 2.5 K and ΔE = 13.3 ± 2.1 kcal/mol for betanova; Tm = 305.6 ± 1.7 K and ΔE = 13.4 ± 1.0 kcal/mol for LLM; Tm = 295.7 ± 3.1 K and ΔE = 9.7 ± 0.5 kcal/mol for beta3s). (b) Nativeness qhb against temperature. Two-state fits were not possible.

These Ehp-based Tm values are high compared to the experimentally detern class="Chemical">mined folded populations, especially for betanova. Note that betanova has a very low hydrophobicity. The correlation between Ehp and folding status is therefore likely to be weak for this peptide.

In contrast to the Ehp-based folded populations, those based on qhb agree quite well with experimental data. In tn class="Chemical">his respect, the situation is the opposite to what we found for the β-hairpins studied above. A possible reason for this difference is discussed below.

3.7 AB zipper

The AB zipper is a designed heterodimeric n class="Chemical">leucine zipper, composed of an acidic A chain and a basic B chain, each with 30 residues [57]. The dimer structure has been characterized by NMR, and a melting temperature of ~340 K was estimated by CD measurements (at neutral pH) [57].

The lowest energy state seen in our simulations is a conformation in which pRMSD calculated over backbone atoms of all ren class="Chemical">sidues in both chains is ~2.7 Å. In this structure, the bRMSD (all residues) of the individual chains A and B to their counterparts in the PDB structure are ~2.5 Å and ~2.4 Å, respectively. Unlike for the other systems described in this article, the boundary conditions have a non-trivial role for this dimeric system. A proper discussion of periodicity, concentration and temperature dependence of this system is beyond the scope of this article. In Fig. 8a, we show the energy landscape, i.e., the mean energy as a function of two order parameters for this system. The X-axis shows the measure pRMSD described earlier. The Y-axis represents the sum of the backbone RMSD of the individual chains. pRMSD can be very large even if the sum of bRMSDs is small: the two chains can be folded without making the proper interchain contacts. Indeed, the figure shows that the major energy gradients are along the Y-axis, showing that it is energetically favorable for both chains to fold to their respective helical states. The correct dimeric native state is energetically more favorable by ~20 kcal/mol compared to two folded helices without proper interchain contacts. This is seen more clearly in Fig. 8b, where we plot the average energy as a function of pRMSD for states with two folded chains. We also simulated the two chains A and B of the dimer in isolation. Both chains folded to their native helical conformations. The melting temperatures estimated based on helix content for chains A and B are 314 K and 313 K, respectively. As indicated above, for the dimer, thermodynamic parameters like Tm cannot be directly estimated from the present simulations.

Figure 8

The heterodimeric AB zipper. (a) Mean energy as a function of pRMSD over both chains and the sum of individual bRMSDs. The direction of the energy gradients implies that a system with two folded monomers is energetically favorable compared to unfolded monomers. The proper dimeric form is the area closest to the origin, and has a lower energy. (b) Mean energy of all states in which both chains have bRMSD < 5 Å, shown as a function of the dimer RMSD measure pRMSD.

3.8 Top7-CFr

Top7-CFr, the C-terminal fragment of the den class="Chemical">signed 93-residue α/β-protein Top7 [58], is the most complex of all molecules studied here. It has both α-helix and β-strand secondary structure elements, and highly non-local hydrogen bonds between the N- and C-terminal strands. CFr is known to form extremely stable homodimers, which retain their secondary structure till very high temperatures like 371 K and high concentrations of denaturants [59].

In [13,14], an earlier version of our model was used to study the folding of CFr. The n class="Chemical">simulations pointed to an unexpected folding mechanism. The N-terminal strand initially folds as a non-native continuation of the adjoining α-helix. After the other secondary structure elements form and diffuse to an approximately correct tertiary organization, the non-native extension of the helix unfolds and frees the N-terminal residues. These residues then attach to an existing β-hairpin to complete the three-stranded β-sheet of the native structure. Premature fastening of the chain ends in β-sheet contacts puts the molecule in a deep local energy minimum, in which the folding and proper arrangement of the other secondary structure elements is hampered by large steric barriers. The above "caching" mechanism, spontaneously emerging in the simulations, accelerates folding by helping the molecule avoid such local minima.

The folding properties of CFr, including the above mentioned caching mechanism, are pren class="Chemical">served under the current modifications of the interaction potential. The centre of the native free-energy minimum shifts from bRMSD (all residues) of 1.7 Å as reported in [13] to about 2.2 Å. This state remains the minimum energy state, although the new energy function changes the energy ordering of the other low energy states. The runs made for this study (see Table 6) found 22 independent folding events. The free-energy landscape observed in the simulations is rather complex with a plethora of deep local minima sharing one or more secondary structure elements with the native structure. They differ in the registry and ordering of strands and the length of the helix. Longer runs are required for the MC simulations to correctly weight these different minima. Temperature dependence of the properties of CFr can therefore not be reliably obtained from these runs.

We note that the simulations ran on twice as many n class="Chemical">processors but were only about one sixth the length of those used for [13], in which 15 independent folding events were found. The improved efficiency is partly due to the changes in the energy function presented here, and partly due to the optimization of the parallel tempering described in [26].

3.9 GS-α3 W

GS-α3 W is a designed n class="Chemical">three-helix-bundle protein with 67 residues [60], whose structure was characterized by NMR [61]. The stability was estimated to be 4.6 kcal/mol in aqueous solution at 298 K, based on CD data [60].

It turns out that this n class="Chemical">protein is very easy to fold with our model. Our results are based on extensive sampling of the conformation space with 64 × 3.5 × 109 Monte Carlo updates, resulting in about 800 independent folding events to the native state. For this estimate, structures with bRMSD (all residues) under 5 Å were taken to be in the native minimum (see Fig. 9 for justification). Two visits to the native state were considered statistically independent (i) if they occurred in independent Markov chains, or (ii) if the two visits to the native state were separated by at least one visit to the highest temperature in the simulation. For the entire run, we spent about 10 days of computing time on 64 AMD Opteron processors running at 2.0 GHz.

Figure 9

The three-helix-bundle protein GS-. (a) Variation of histogram of bRMSD with temperature. At high temperatures, there is a broad distribution of bRMSD with values > 10 Å. At lower temperatures there are three clearly separated clusters. Representative structures from these clusters are also shown (color) aligned with the native structure (gray). (b) Temperature dependence of specific heat, C, and the ratio hof the observed helix content and the helix content of the native structure.

The three-helix-bundle n class="Chemical">protein GS-. (a) Variation of histogram of bRMSD with temperature. At high temperatures, there is a broad distribution of bRMSD with values > 10 Å. At lower temperatures there are three clearly separated clusters. Representative structures from these clusters are also shown (color) aligned with the native structure (gray). (b) Temperature dependence of specific heat, C, and the ratio hof the observed helix content and the helix content of the native structure.

In Fig. 9a, we show how the probabilities for structures with different bRMSD vary with temperature in the n class="Chemical">simulations. Clearly, the protein makes a transition from a rather continuous distribution of bRMSD at high temperatures to a distribution dominated by three well separated clusters. Analysis of the structures at the lower temperatures shows that all three free-energy minima consist almost exclusively of structures with all three helices of GS-α3 W formed. The plot of the ratio of the observed helix content and the helix content of the native state, shown in Fig. 9b, further supports this idea. The average value of this ratio approaches 1 as the temperature decreases below 300 K. The specific heat curve, also shown in Fig. 9b, indicates that the formation of these structures correlates with the steepest change in energy.

The cluster with a center at bRMSD ~3 Å dominates at the lowest temperatures. The structures contributing to the cluster with ~8–9 Å bRMSD superficially look like well folded n class="Chemical">three-helix bundles. But as illustrated in the figure, the arrangement of the helices is topologically distinct from the native arrangement. The cluster seen at larger bRMSD values is broader and consists of a host of structures in which two of the helices make a helical hairpin, but the third helix is not bound to it. The unbound helix could be at either side of the chain.

According to our model therefore, the population at the lowest temperatures consists of ~80% genuinely native structures, ~10% n class="Chemical">three-helix bundles with wrong topology, and ~10% other structures with as much helix content as the native state. In order to experimentally determine the true folded population of the protein, the experimental probe must be able to distinguish the native fold from the other helix rich structures described here.

4 Discussion

The model presented here is intrinsically fast compared to many other all-atom models, because all interactions are short range. By exploiting tn class="Chemical">his property and using efficient MC techniques, it is possible to achieve a high sampling efficiency. We could, for example, generate more than 800 independent folding events for the 67-residue GS-α3 W. The speed of the simulations thus permits statistically accurate studies of the global free-energy landscape of peptides and small proteins.

In developing this potential, a set of 17 n class="Chemical">peptides with 10–37 residues was studied. The peptides were added to this set one at a time. To fold a new sequence sometimes required fine-tuning of the potential, sometimes not. A change was accepted only after testing the new potential on all previous sequences in the set. In its final form, the model folds all 17 sequences to structures similar to their experimental structures, for one and the same choice of potential parameters.

Also important is the stability of the peptides. A small polypeptide chain is unlikely to be a clear two-state folder, and therefore its apparent folded population will generally depend on the obn class="Chemical">servable studied. For β-sheet peptides, we used the hydrophobicity energy Ehp and the hydrogen bond-based nativeness measure qhb to monitor the melting behavior. The extracted Tm values indeed showed a clear probe dependence; the Ehp-based value was always larger than that based on qhb. For the β-hairpins studied, we found a good overall agreement between our Ehp-based results and experimental data. For the three-stranded β-sheets, instead, the qhb results agreed best with experimental data. The reason for this difference is unclear. One contributing factor could be that interactions between aromatic residues play a more important role for the β-hairpins studied here than for the three-stranded β-sheets. These interactions may influence spectroscopic signals and are part of Ehp. Probe-dependent Tm values have also been obtained experimentally, for example, for trpzip2 [53].

The probe dependence makes the comparison with experimental data less straightforward. n class="Chemical">Nevertheless, the results presented clearly show that the model captures many experimentally observed stability differences. In particular, among related peptides, the calculated order of increasing thermal stability generally agrees with the experimental order, independent of which of our observables we use.

It is encouraging that the model is able to fold these 17 sequences. However, there is no existing model that will fold all peptides, and our model is no exception. Two sequences that we unsuccessfully tried to fold are the β-hairpins n class="Chemical">trpzip4 and U16, both with 16 residues. Trpzip4 is a triple mutant of GB1p with four tryptophans [52]. For trpzip4, our minimum energy state actually corresponded to the NMR-derived native state [52], but the population of this state remained low at the lowest temperature studied (~14% at 279 K, as opposed to an estimated Tm of 343 K in experiments [52]). U16 is derived from the N-terminal β-hairpin of ubiquitin [62]. It has a shortened turn and has been found to form a β-hairpin with non-native registry [62]. In our simulations, this state was only weakly populated (~8% at 279 K, as opposed to an estimated ~80% at 288 K [62]). Instead, the main free-energy minima corresponded to the two β-hairpin states with the registry of native ubiquitin, one with native hydrogen bonds and the other with the complementary set of hydrogen bonds.

Our calibration of the potential relies on experimental data with non-negligible uncertainties, on a limited number of n class="Chemical">peptides. It is not evident that this potential will be useful for larger polypeptide chains. Therefore, as a proof-of-principle test, we also studied three larger systems, with very good results. Our simulations showed that, without having to adjust any parameter, the model folds these sequences to structures consistent with experimental data. Having verified this, it would be interesting to use the model to investigate the mechanisms by which these systems self-assemble, but such an analysis is beyond the scope of this article. The main purpose of our present study of these systems was to demonstrate the viability of our calibration approach.

The potential can be further constrained by confronting it with more accurate experimental data and data on new sequences. The challenge in this n class="Chemical">process is to ensure backward compatibility – new constraints should be met without sacrificing properties already achieved.

5 Conclusion

We have described and tested an implicit solvent all-atom model for n class="Chemical">protein simulations. The model is computationally fast and yet able to capture structural and thermodynamic properties of a diverse set of sequences. Its computational efficiency greatly facilitates the study of folding and aggregation problems that require exploration of the full free-energy landscape. A program package, called PROFASI [28], for single- and multi-chain simulations with this model is freely available to academic users.