Folding free energy landscapes of β-sheets with non-polarizable and polarizable CHARMM force fields

The so-called “protein folding problem,” namely, how proteins arrive at a particular three-dimensional shape determined primarily by their sequence of amino acids, remains one of the grand challenges of molecular biology.¹ While most estimates of the number of distinct protein folds are around 10³ (see Ref. 2), they are practically all derived from the two most common secondary-structure elements: the α-helix and the β-sheet. Thus, model systems exemplifying one of these elements represent important tools for the study of protein folding.

One of the most common representatives of the β-sheet used in protein folding studies is the B1 domain of streptococcal protein G (GB1) as well as its C-terminal β-hairpin (residues 41-56), which comprises just 16 amino acids. Although other isolated β-hairpin peptides have been found,^3,4 GB1 has proven particularly amenable to experiments looking at folding dynamics and mechanisms.⁵ GB1 is stabilized by six native backbone hydrogen bonds and a hydrophobic core consisting of four residues (see Fig. 1). Muñoz et al. measured folding to be a two-state process taking ∼6 μs, and multiple studies have found the fraction folded to be 30%–50% at room temperature.^5,6 This folding time is ∼30× slower than that for α-helices.^7,8 The earliest folding model put forward proposed that the β-turn forms first and then the hairpin “zips up” from there, with the hydrophobic core forming last.^5,9

Because of its small size and fundamental nature, the folding of the GB1 β-hairpin has inspired a number of molecular dynamics (MD) studies. In particular, multiple groups have calculated the free-energy landscape of GB1. For example, Dinner et al. chose two reaction coordinates to map GB1’s conformations, namely, the native hydrogen bonds and the packing of the hydrophobic core residues; they found, contrary to the earlier model, that the hydrophobic core forms first, followed by hydrogen bond formation along the backbone.¹⁰ Successive simulation studies agreed with this model.^11,12 Zhou and Berne compared explicit and implicit solvent models with the same force field, finding that the folded probability of GB1 in the implicit solvent was much lower than that in the explicit solvent, the latter matching experiments.¹³ More recently, Shao et al. demonstrated that a particular implicit solvent model, GB^OBC, when used with the AMBER ff96 force field, could reproduce the expected structures of GB1.¹⁴ Thus, GB1 β-hairpin folding developed into a useful and sensitive test of simulation force fields. Best and Mittal further demonstrated this utility by comparing six force-field/explicit-water-model combinations; they found that the primary differences between force fields are the residual secondary structures observed in the unfolded state as well as the dimensions of the unfolded state, which was often too compact.¹⁵ The delicate balance between solute-solvent and solute-solute interactions was cited as a reason for the differences.

Force fields that include explicit polarizability have grown in popularity over the years, in part due to improvements in the methods and also the expansion in computational resources available.^16–19 Two commonly used approaches that have been extended to proteins are the AMOEBA force field²⁰ and the Drude-2013 force field.²¹ The AMOEBA-2013 force field proved adept at maintaining conformations of ten small proteins with varied folds (α-helices, β-sheets, and mixed) over 30-ns simulations.²⁰ The Drude-2013 force field was used for simulations of a few small peptides as well as ten proteins, finding them to be stable on the 100-ns time scale.²¹ Dynamic polarization of the peptide backbone led to improvements in both the temperature dependence and cooperativity of α-helix folding for the (AAQAA)₃ peptide.²² Lopes et al. also simulated the GB1 β-hairpin with Drude-2013. While it remained close to the crystal structure over the course of a 150-ns simulation, the comparison of estimated NMR chemical shifts with experiment for residues 42-45 proved worse than the non-polarizable CHARMM36 force field.²¹ Additionally, backbone hydrogen bonds in β-sheet structures were approximately 0.05 Å too long compared to high-resolution crystal structures.²¹ 

In the present study, we determine the free energy of GB1’s conformational landscape for three force fields: CHARMM36 (Ref. 23), CHARMM22^* (Ref. 24), and the 2015 release of the Drude-2013 model (Ref. 21). We find that while the β-hairpin of GB1 is stable in the non-polarizable models, it is unstable in the Drude-2013 model, with the latter favoring an unfolded, extended state. Hydration free energies reveal that amide-containing groups are more soluble in the Drude-2013 model than in the non-polarizable models, resulting in a more soluble backbone. However, this enhanced solubility, which is still less than that found in experiments, does not appear to be countered by a corresponding enhancement in peptide–peptide interactions, namely, for the native hydrogen bonds that form the β-hairpin. Using a perturbative approach, we show that small increases to the backbone amide N polarizability can improve both the backbone solubility and folding properties of the GB1 β-hairpin in the Drude-2013 model without sacrificing the improvements observed for α-helical peptides.

All-atom simulations of non-polarizable GB1 were performed starting from residues 41–56 of PDB 1GB1 (Ref. 25; see Fig. 1), capped with an acetylated N-terminus and amidated C-terminus. For non-polarizable models, the visualization and analysis program VMD²⁶ was used to place GB1 in a water box of 14 314 TIP3P (Ref. 27) water molecules with dimensions 66 × 66 × 104 ų. Systems were neutralized with 6 sodium and 3 chloride ions, as used previously,¹⁵ for a total of 43 204 atoms. A harmonic restraint was added to keep the end-to-end vector along the z-axis. Molecular dynamics simulations were carried out using NAMD 2.10-12 (Ref. 28) with CHARMM all-atom force fields C22^* (Ref. 24) and C36 (Refs. 23 and 29).

The temperature was fixed at 300 K using Langevin dynamics; the pressure was kept constant at 1 atm using the Langevin piston method.³⁰ The equations of motion were integrated using the RESPA multiple time step algorithm with a time step of 2 fs used for all bonded interactions, 2 fs for short-range non-bonded interactions and 4 fs for long-range electrostatic interactions. Long-range electrostatic interactions were calculated using the particle-mesh Ewald method with a real-space cutoff of 12 Å.³¹ Short-range non-bonded Lennard-Jones interactions were cutoff at 12 Å with a potential switching function beginning at 10 Å bringing the potential energy smoothly to zero at the cutoff distance. Bonds involving hydrogen atoms were constrained to their equilibrium length, employing the SETTLE algorithm³² for water molecules and the SHAKE algorithm for all others.³³ 

For the polarizable model, a non-polarizable model was first built using CHARMM-GUI.³⁴ Similar to the previous non-polarizable models, GB1 was solvated in a water box of dimensions 66 × 66 × 104 ų with 14 256 water molecules and neutralized with 6 sodium and 3 chloride ions. An additional harmonic restraint was added to keep the end-to-end vector along the z-axis. The Drude Prepper³⁴ from CHARMM-GUI was then used to convert the non-polarizable model into the Drude-2013 polarizable model for the protein²¹ and ions³⁵ solvated in SWM4-NDP polarizable water molecules,³⁶ resulting in a system of 71 748 atoms, including Drude particles. We used the Drude-2013 parameters released in July 2015. The system was then minimized and pre-equilibrated using the NAMD input scripts provided by Drude Prepper.

The temperature for parent atoms was fixed at 300 K using Langevin dynamics; a separate Langevin thermostat was coupled to the Drude particles with temperature 1 K. Additionally, a hard wall restraint at 0.25 Å was added for all Drude particle-parent atom bond lengths. Thole corrections to electrostatic interactions³⁷ were also extended to nonbonded pairs of Drude oscillators for which a Thole screening parameter is defined within a 5.0-Å cutoff. Analytic long-range Lennard-Jones corrections were also implemented.³⁸ The pressure was kept constant at 1 atm using the Langevin piston method.³⁰ The equations of motion were integrated using a 1-fs time step for all interactions. Long-range electrostatic interactions were calculated using the particle-mesh Ewald method;³¹ cutoffs were the same as for the non-polarizable simulations above. Bonds involving hydrogen atoms were constrained only for water molecules.

Two-dimensional potentials of mean force (PMFs) were calculated using umbrella sampling with replica-exchange (REUS).³⁹ Reaction coordinates were calculated and biased using the collective variable (colvars) module of NAMD 2.10 (Ref. 40). The first reaction coordinate is the number of native hydrogen bonds (N_hb), calculated using the hbonds colvar. N_hb is the sum of individual dimensionless hydrogen bond scoring functions, hbf, calculated for each donor and acceptor pair,

where x_i is the Cartesian coordinates of atom i and d₀ = 3.3 Å is the distance cutoff for an occupied hydrogen bond. hbf ∈ [0, 1], with values near 0 for donors and acceptors far outside the cutoff distance and values near 1 for those well inside the cutoff. Therefore, N_hb ∈ [0, 6] for the six native hydrogen bonds of GB1. The second reaction coordinate is the radius of gyration of the hydrophobic core (R_G), i.e., the non-hydrogen atoms of residues Trp43, Tyr45, Phe52, and Val54. Either 89 (C36 and Drude-2013) or 100 (C22^*) windows were simulated for 15 ns/window, for a total simulation time of 1.34-1.50 μs per system. Windows were spaced along the N_hb coordinate by 0.25 for N_hb ≤ 1.5 and by 0.5 for N_hb > 1.5, using a 25 kcal/mol harmonic force constant. Windows were spaced along the R_G coordinate by 1.0 Å, using a 6.25 kcal/mol·Å² harmonic force constant. The first 2 ns/window of each system was omitted before calculating the PMF, generated via the weighted histogram analysis method (WHAM).^41,42 Starting states for each window were generated using steered molecular dynamics (SMD).⁴³ 

1D PMFs were calculated by integrating out the second coordinate as follows:⁴⁴ 

where W(x, y) is the 2D PMF, $w$(x) is the corresponding 1D PMF, and (x_c, y_c) is an arbitrary point in the collective variable space. To calculate folding free energies, ΔG_fold, from the 1D PMFs, we first located the peak of the PMF between the two extrema of the x-coordinate; then, for folding-related coordinates (e.g., the number of hydrogen bonds or the fraction of the native structure), the 1D PMF was integrated to the left and to the right of the peak to get G_unfolded and G_folded, respectively, where ΔG_fold = G_folded − G_unfolded. For distance-related coordinates (e.g., the radius of gyration or end-to-end distance), the definitions for folded and unfolded are reversed. Errors for ΔG_fold were calculated as the standard deviation of ΔG_folds calculated from 1-ns/window blocks of the REUS trajectories.

Drude particle charges, q_D, and parent atom charges, q_P, are determined by

where q_A is the atomic partial charge, K_D = 500 kcal/mol·Å² is the strength of the harmonic potential connecting the Drude particle to its parent atom, α is the atomic polarizability, and sgn(α) = +1 if α > 0, −1 if α < 0, and 0 if α = 0. Parent charges also include any associated lone pair particles. If we modify the polarizability by a factor α′/α = γ, then $qD′$ = γ^1/2q_D and $qP′$ = q_A − $qD′$⁠. For atoms with lone pair particles, we keep the parent charge constant and split the new charge between the lone pair particles.

States generated from the REUS simulations were used to generate new PMFs perturbatively. First, the total potential energy was calculated for 10-ps snapshots from the last 5 ns/window of the REUS simulations. Then, new potential energies were calculated for these same snapshots using adjusted parameters. WHAM histograms were reweighted by the difference in the energies using the following equation, and new PMFs were calculated:

where U_new and U_old refer to the potential energy U(r(t)) with the new and old parameters, respectively; r(t) is the value of the collective variables at time step t; and r_ij is the value of the collective variables in bin ij. For consistency, the original PMFs were also recalculated using the same 10-ps snapshots to account for any bias in the reduced data set.

The sidechain hydration energies were calculated using CHARMM c40b1.⁴⁵ The 40-stage Weeks–Chandler–Anderson (WCA)-decomposition free-energy perturbation (FEP) procedure described by Deng and Roux was used.^46,47 The electrostatic and dispersion components were calculated from windows from λ = 0.0, 0.1, 0.2, …, 1.0. Exchanges were performed between neighboring replicas with a frequency of 2 ps. The solutes were placed in a 32 × 32 × 32 Å simulation cell of TIP3P-model water (roughly 980 water molecules). The Gibbs energy calculations were performed using a 2 ns simulation for equilibration followed by a 5 ns simulation to sample the free energies. Gibbs energies were calculated using WHAM.⁴¹ For charged amino acids, the electrostatic component solvation energies were corrected for an interfacial potential of −520 mV.⁴⁸ Uncertainties were estimated by calculating the standard deviation of the free energies taken from thirds of the production simulation.

A correction for the neglect of dispersion interactions outside the 12-Å cutoff distance was determined by calculating the average of the difference between the interaction energy with a 12-Å Lennard-Jones cutoff and that with a 45-Å Lennard-Jones cutoff in a 1-ns MD simulation. This energy is included in the dispersion component of the hydration energy.

The Drude-2013 model hydration energies were calculated using the same procedure; however, a 1-fs time step was used. The SWM4-NDP water model was used as the solvent. The electrostatic components of the FEP of the Drude-2013 hydration energies were determined without exchanges. An interfacial potential of −545 mV was used to calculate the correction to the electrostatic component of the solvation energy for ionic solutes.⁴⁹ 

Two-dimensional potentials of mean force (PMF) were calculated for GB1 using CHARMM36 (C36), CHARMM22^* (C22^*), and CHARMM Drude polarizable (Drude-2013) force fields (see Fig. 2) using replica exchange umbrella sampling (REUS; see Sec. II). Previous simulation studies of GB1 near room temperature using a variety of protein force fields and solvent models have found a number of local free energy minima, including the fully folded states, partially folded states, compact unfolded states, and extended states.^{11,13–15,50–54} We find that the three CHARMM force fields also produce a variety of local minima, resulting in distinct folding energetics for each model. The free energy landscape for C36 shows the native β-hairpin (∼4–5 native hydrogen bonds) to be very stable, roughly 6 kcal/mol lower in energy than the compact, unfolded state and 8 kcal/mol lower than the extended state. The folded state is much less stable in C22^*, with a free energy basin composed of multiple local minima extending from the native hairpin with 5 hydrogen bonds to a partially unfolded hairpin with 2 hydrogen bonds, with the lowest energy state residing near 3 hydrogen bonds; local minima are all separated by barriers of <2 kcal/mol. This free energy basin is around 2–4 kcal/mol lower in energy than both the compact and extended unfolded states. While C36 and C22^* have global minima in the native and partially folded states at coordinates (N_hb, R_G) = (4.95, 5.45) and (3.11, 5.45), respectively, the Drude-2013 model shows a global minimum in the extended, unfolded state at coordinates (0.23, 11.65), with shallower local minima in the partially folded and fully folded regions.

To calculate the folding energetics, we integrated out the N_hb and R_G coordinates separately from our 2D PMFs to generate 1D PMFs along the other corresponding coordinate (see Fig. S1 of the supplementary material) and calculated the free energy of folding, i.e., ΔG_fold = G_folded − G_unfolded (see Sec. II). The results are shown in Table I. From the 2D PMFs, the N_hb coordinate appears to better differentiate between the folded and unfolded states than the R_G coordinate, with a clear barrier between the two states for all three force fields. Hereafter, we will only compare folding free energies along the N_hb coordinate. GB1 is between 30% and 50% folded at room temperature,^5,6 so ΔG_fold = 0.0–0.5 kcal/mol. C36 overstabilizes the folded state by ∼8 kcal/mol, C22^* overstabilizes it by ∼2 kcal/mol, and Drude-2013 destabilizes the folded state by ∼2 kcal/mol. ΔG_fold is within the error of the experimental range for C22^* and is just outside this range for the Drude-2013 model, but is far too negative for C36. Therefore, the C22^* and Drude-2013 models are roughly equivalent in describing the folding behavior of the GB1 β-hairpin, with the C36 model performing the worst of the three force fields.

To examine the differences in the GB1 stability between the three force fields, we calculated hydration energies of amino-acid side-chain analogs using the C36 and Drude-2013 models (Table II). With the exception of Asp, Glu, Arg, and Trp side-chain partial charges, C22^* non-bonded parameters are equivalent to C36.^24,55 In the analogs used, the side chain is terminated with a hydrogen atom at the position where the α-carbon would be bonded. Molecules constructed to mimic the side chain of an amino acid do not necessarily serve as an appropriate measure of the hydration of a side chain attached to a protein,⁵⁶ but they do serve to test if the non-bonded parameters for the chemical moieties present in the side chains are appropriate for describing their interactions with water.^57–59 

In general, hydration energies predicted by both the C36 and Drude-2013 models are in good agreement with the experimental hydration energies. The C36 and Drude-2013 models both show a significant improvement in Trp side chain hydration over C22 (equivalent to C22^*), which underestimates the free energy by 3 kcal/mol,⁴⁶ whereas the C36 and Drude-2013 models only underestimate it by 0.5 and 0.8 kcal/mol, respectively. The Drude-2013 model also shows a systematic improvement for sulfur-containing side chains, Cys and Met, over the C36 model. The Drude-2013 model correctly predicts these to be sparingly soluble, while the C36 model underestimates their solubility by 1 and 2 kcal/mol, respectively. Conversely, the hydration energies of non-aromatic hydroxyl containing groups (i.e., Ser and Thr) are underestimated by the Drude-2013 model, while the C36 model predicts them accurately. This difference appears to stem from the electrostatic component of the hydration energy: in the Drude-2013 model, they are roughly 1 kcal/mol smaller than those of the C36 model (see Tables SI and SII of the supplementary material). This is consistent with the O—H bond in the Drude-2013 model being less polar (q_H = 0.36e in the Drude-2013 model, while q_H = 0.43e in the C36 model). Polar and non-polar residues alternate along the two strands of the GB1 β-hairpin, a common motif in anti-parallel β-strands,⁶⁰ which situates the polar and non-polar residues to point in the opposite direction, leading to distinct polar and non-polar faces of the β-sheet. Given that four of the six hydrogen-bonding residues in GB1 are threonines, reduced solvation of the threonine side chains in the Drude-2013 model could have a compounding effect, reducing the hydration free energy of the polar face of the β-hairpin.

The C36 model underestimates the hydration energies of amides (i.e., NMA, Gln, and Asn) by roughly 2 kcal/mol. One justification for this underestimation is that the force field compensates for the lack of induced polarization in the model. When a protein folds into secondary structures, hydrogen bonding between backbone amides results in polarization of the C=O and N—H bonds, which gives a complementary effect that favors the propagation of secondary structures.²² Non-polarizable force fields like C36 are incapable of describing this effect rigorously, but by defining the water–amide Lennard-Jones parameters such that the interactions are weakened, the model has the correct propensity for forming secondary structures.

In contrast to the C36 model, the Drude-2013 model hydration energies of the amide-containing molecules (i.e., NMA, Gln, and Asn) are in better agreement with the experimental values. This suggests a possible explanation for the Drude-2013 model predicting a (slightly) unstable β-hairpin, as opposed to a stable β-hairpin for C36. The amide hydration energies of the Drude-2013 model are more negative than the C36 model by roughly 1 kcal/mol. If the stabilization that is achieved by the backbone–backbone hydrogen bonds of the β-hairpin does not compensate for the increased hydration energy of the amides, the unfolded state will be favored.

We observed key differences in hydration free energies between the Drude-2013 and C36 models that lead to distinct free energy profiles for GB1 β-hairpin folding. The Drude-2013 model improves the hydration free energies over C36 for some chemical motifs, namely, amides and sulfur-containing compounds, while it underestimates the free energies for other motifs, such as hydroxyl groups. Although the NMA hydration free energy is improved in the Drude-2013 model, it is still underestimated by ∼1.5 kcal/mol relative to the experimental value. Enhancing the atomic polarizability of backbone amide N and side chain hydroxyl O atoms should improve the hydration of the peptide backbone and Thr and Ser side chains, respectively. We recalculated the electrostatic component of the hydration free energies for NMA, Thr, and Ser for atomic polarizabilities of amide N and hydroxyl O atoms scaled by up to 2× (see Fig. 3).

For NMA, experimental hydration free energies are most accurately reproduced when the polarizability of the amide N atom is increased by 60%. For Ser and Thr side chains, increasing the polarizability of the hydroxyl O atom by 50% and 30%, respectively, was the most accurate in reproducing experimental hydration free energies. For all of NMA, Ser, and Thr, increasing their respective atomic polarizabilities by 40% reproduces experimental hydration free energies within ∼0.5 kcal/mol. Adjusting the atomic charge of the amide H atom to its C36 value also reproduced NMA hydration free energies fairly well, underestimating the free energy by less than 0.5 kcal/mol; adjusting the amide N atomic charge to its C36 value, however, overestimates the hydration free energy by more than 1 kcal/mol (see the supplementary material). Overall, side chain hydroxyl free energies are more sensitive to changes in polarizability and dipole moment than amides. This might be expected since, although N atoms are more polarizable than O atoms, O—H groups have larger intrinsic dipole moments than N—H groups. In addition, C=O groups have a large effect on the solubility of amides because they have larger dipole moments, both intrinsically and when hydrogen bonded to water, than the N—H groups (see Fig. S6 of the supplementary material).

We recalculated our folding PMFs for GB1 using these new polarizability parameters (see Sec. II and Fig. 4). We extended the range of polarizabilities tested for the N and O atoms down to 0 in order to better understand the role of polarizability in the folding behavior. As might be expected, the folded state becomes more favorable as the backbone amide N atom becomes more polarizable, with ΔG_fold decreasing roughly monotonically with increasing polarizability, becoming negative when the polarizability is increased by 70% (1.7×) compared to the original Drude-2013 parameter. For increases of 30%–60% (1.3–1.6×), ΔG_fold decreases by ∼1.4–1.9 kcal/mol over the original polarizability, with 30% showing the largest drop, bringing ΔG_fold closer to the experimental values of 0.0–0.5 kcal/mol.^5,6 A similar shift is seen when adjusting the atomic partial charges of the amide N and H atoms to their C36 values (see the supplementary material). Because using reduced data sets can introduce a sampling bias to the free energy (see Sec. II), we approximate this bias by subtracting the full-data free energy from the reduced-data free energy: $ΔGfoldbias=ΔGfoldreduced−ΔGfoldfull$ = +1.4 kcal/mol. If we adjust for this bias, increasing the polarizability 30%–60% should produce ΔG_fold values only ∼0.5–1.0 kcal/mol greater than experimental values, which is within the error of the PMF calculations (see Table I). However, the change in parameters could be affected differently by the sampling bias, meaning that full REUS calculations would be needed to verify the improvement in the free energies. From our perturbative approach, either a 30%–60% increase in the amide N polarizability or a reversion of the backbone amide H atomic charge to its C36 value appears to be a viable means of improving both the backbone solubility and β-hairpin folding properties in the Drude-2013 model. Alternatively, decreasing the Lennard-Jones minimum between the backbone amide N and carbonyl O atoms by 7% also reproduces experimental folding free energies, but without the additional improvement in solubility (see the supplementary material).

Unlike the backbone amide N polarizabilities, altering side chain hydroxyl O polarizabilities does not monotonically change the folding free energy of GB1. GB1 contains five Thr residues, four of which are involved in native β-hairpin hydrogen bonding with the fifth residing in the turn region. Increasing the polarizability of Thr side chain hydroxyl O atoms by 30%, which most accurately reproduces experimental hydration free energies [see Fig. 3(b)], results in the largest, positive folding free energy for all tested polarizabilities, shifting it to the unfolded state by +4.4 kcal/mol over the original Drude-2013 parameters. At larger polarizabilities, the folding free energy stabilizes to around 1 kcal/mol lower than the original polarizability. On the other hand, reducing the polarizability of the hydroxyl O atom close to 0, i.e., making the Thr side chains more hydrophobic, stabilizes the folded state to the same extent as greatly increasing the backbone amide N polarizability. The increase in folding free energy around 1.3× the original O polarizability appears to be driven by an increase in the free energy of the intermediate region between N_hb = 2–4. For higher polarizabilities, this is further accompanied by a corresponding decrease in the free energy near the folded state. A more modest increase of 20% in the polarizability shifts the folding free energy up by 1 kcal/mol over the original parameters while underestimating the Thr side chain hydration free energy by only 1 kcal/mol. This appears to be a reasonable trade off, especially if combined with the aforementioned adjustments to the backbone amide polarizability or atomic partial charges.

We tested the proposed changes to the Drude-2013 parameters for an α-helical peptide, Ala₁₀, a decamer of alanines (see the supplementary material for more details). Ala₁₀ folding free energy landscapes for the three CHARMM force fields tested here follow roughly the same patterns observed for GB1 folding, namely, that the Drude-2013 model favors the unfolded, extended state more so than the additive C36 and C22^* models, with ΔG_fold > 0 for all three models (see Fig. S8 and Table S3 of the supplementary material). For the α-helical peptide, (AAQAA)₃, the Drude-2013 model correctly predicted a positive folding free energy (only ∼25% folded), although this was actually slightly higher than for C36 (∼20% folded).²² For Ala₁₀, we can roughly compare the folding free energies to those of Ala₅. Graf et al. showed that backbone NMR coupling parameters were roughly constant for Ala₃ up to Ala₇, and fitting MD simulations to the NMR results showed very little change in the distribution of structures between Ala₅ and Ala₇, with no α-helical content in either.⁶³ Best et al. repeated this analysis for Ala₅ with a larger sample of MD simulations using multiple force fields and found it to contain ∼2%–11% helical content (ΔG_fold ∼ 1.2–2.3 kcal/mol). For Ala₁₀, the Drude-2013 model sits right at the upper bound of this range, while C36 and C22^* sit just below this range, although they are all within one standard deviation of it. Only small increases in amide N polarizabilities of 10%–20% marginally decrease ΔG_fold, while all others cause ΔG_fold to increase. We see similar results for the Lennard-Jones parameters, where a decrease in the N⋯O LJ minima of up to 4% yields small decreases in ΔG_fold, with all larger deviations increasing ΔG_fold. C36 backbone amide charges also leave the folding free energy outside of the Ala₅ range. This suggests, along with good experimental agreement for (AAQAA)₃, that the Drude-2013 parameters are optimized for α-helices, but may be less so for β-sheets. A modest increase in backbone amide N polarizabilities of 30% could be sufficient for improving both the backbone hydration and β-sheet folding free energies, while not appreciably affecting α-helical folding free energies.

Separate from the issues with folding GB1, the amino-acid-analog hydration energies suggest other areas where the parameters for the Drude-2013 force field could be improved. As noted above, the hydration energies of the hydroxyl-containing amino acids are underestimated by the Drude-2013 model. The parameters for the charged side chains could also be improved; the hydration energies of Glu and Asp are overestimated, while the hydration energy of Lys is underestimated. Although there is no experimental value to compare to, the Drude-2013 model hydration energy for Arg is 4 kcal/mol lower than predicted by the C36 model.

The Drude-2013 model hydration energy of Asn is more negative than that of Gln (−9.7 kcal/mol vs. −8.1 kcal/mol), although the experimental hydration energies are essentially equal (−9.7 kcal/mol and −9.4 kcal/mol, respectively). This difference results from the use of a smaller Lennard-Jones radius for the C_β of Asn (atom type CD32C, σ = 3.22 Å) in comparison to Gln (atom type CD32A, σ = 3.74 Å); when the hydration energy of Gln is calculated using the CD32C atom type for C_β, it becomes −11.2 kcal/mol. The σ parameter for the Asn CD32C atom type is actually anomalously small for an aliphatic carbon, which results in a smaller repulsive component of the hydration energy. This suggests that the Drude-2013 model hydration energy of Asn only agrees with the experimental value because of the small radius of C_β, but it does not accurately represent the strength of the water–amide interactions.

We also examined interactions in the β-turn region of GB1 and how they correlate to folding. Salt bridges that stabilize the turn have been shown to increase the stability of the β-hairpin for another model β-sheet peptide, the PIN1 WW domain.⁶⁴ For GB1, however, the Asp47–Lys50 salt bridge in the turn region, although appearing in some crystal structures, does not consistently form in solution.⁶⁵ For C22^* and Drude-2013, the salt bridge does not form, while for C36, it forms before hydrogen bonding begins and is maintained throughout the folding process (see Fig. S10A of the supplementary material). This salt bridge appears to be the source of the overstabilized β-hairpin in the C36 model. In the Drude-2013 model, Lys50 dipole moments are significantly reduced compared to C36 (a decrease of ∼3–4 D; see Fig. S5 of the supplementary material), while other non-glutamate residues remain roughly the same between the three force fields (see Fig. S5 of the supplementary material). Interestingly, Asp47 dipole moments are also significantly reduced in C22^* compared to C36, most likely due to decreased partial charges of the aspartate side-chain carboxylate group in C22^* (Ref. 24) as a similar decrease in dipole is also observed for Asp46 (see Fig. S5 of the supplementary material). The decrease in dipole moments for Asp47 and Lys50 in the C22^* and Drude-2013 models, respectively, from their C36 levels could explain why the salt bridge is only observed in the C36 model.

Finally, we wanted to address experimental measures of protein folding and how this compares to folding observed in MD simulations. Muñoz et al. measured GB1 folding using fluorescence measurements of the buried Trp43 side chain.⁵ The solvent accessible surface area (SASA) correlates well with tryptophan fluorescence,⁶⁶ so we calculated the SASA of Trp43 for our REUS simulations. As one might expect, Trp43 SASA correlates well with the R_G coordinate but poorly with the N_hb coordinate (see Fig. 5). If we integrate the 2D PMFs over the fraction of the exposed Trp43 surface area, we see only a single well for all three models (see Fig. S11 of the supplementary material). This follows what was observed for the R_G coordinate, where only the original Drude-2013 model showed two-state behavior (see Fig. S1 of the supplementary material). Essentially, since formation of the hydrogen bonds follows the hydrophobic collapse, the exposed tryptophan surface area only tracks the hydrophobic collapse portion of the folding process, giving no indication of hydrogen bond formation. This agrees well with experiments that have found that folding rates determined by tryptophan fluorescence measurements are often faster than the true full folding rate.⁶⁷ Thus, tryptophan fluorescence alone is not a sufficient reaction coordinate for protein folding and, instead, requires at least one supplemental coordinate to capture the complete folding process.

We calculated folding free energy landscapes of the GB1 β-hairpin using CHARMM non-polarizable force fields and the CHARMM Drude-2013 polarizable model and found that the Drude-2013 model destabilizes the β-hairpin, resulting in ΔG_fold ∼ +2 kcal/mol, versus −7 and −2 kcal/mol for C36 and C22^*, respectively. Hydration free energies show improved solvation of the peptide backbone in the Drude-2013 model over the non-polarizable C36, suggesting that peptide–peptide hydrogen bonding is too weak to compensate for a more solvated backbone. By increasing the backbone amide N atomic polarizabilities by 30%, we were able to further improve the backbone solubility as well as shift the folding free energy closer to the experimental value of ∼+0.5 kcal/mol.^5,6

Improved backbone hydration energies in the Drude-2013 model may be the result of numerous occurrences. First, C36 amide hydration energies are underestimated by ∼2–3 kcal/mol (see Table II). This presumably compensates for the lack of N—H and C=O polarization during hydrogen bonding. However, although this works well for α-helices,⁶⁸ it results in overstabilized β-sheets. Altered dihedral parameters and the removal of the helix-overstabilizing CMAP potential allow C22^* to avoid such a strong preference for the folded state.^24,69 For the Drude-2013 model, altered partial charges combined with dynamic atomic polarization strengthen backbone–water hydrogen bonding. O partial charges are more negative in the SWM4-DNP water model (−0.834e for TIP3P compared to −1.115e for SWM4-NDP), and they have a deeper Lennard-Jones potential well (−0.15 kcal/mol for TIP3P compared to −0.21 kcal/mol for SWM4-NDP), both of which are typical for 4-point water models.^70–72 Additionally, the net charge for backbone N—H groups increases from −0.16e in C36 to −0.11e in Drude-2013, while the net charge for backbone C=O groups decreases from 0e in C36 to −0.04e in Drude-2013. The reduced charge for N—H groups is due to a less negatively charged N atom and a less positively charged H atom, while the net negative charge for C=O groups is a result of a less positively charged C atom and a more negatively charged O atom, with the negative charge on the O atom split between two lone pair charges. N—H and C=O dipole moments are also typically higher in the Drude-2013 model than in the C36 or C22^* models. For native hydrogen-bonding residues, N—H dipole moments in the Drude-2013 model are enhanced by the polarization of the N atoms and alignment of the dipoles during folding. Conversely, C=O dipole moments are enhanced primarily by the alignment of the atomic dipoles, as the C and O atoms are less polarizable than the N atoms. However, C and O atoms do polarize to a larger extent when involved in water–peptide hydrogen bonding. As we also observed with hydroxyl groups in polar side chains, the reduced partial charges on the backbone N and H atoms in the Drude-2013 model may be decreasing its affinity for hydrogen bonding. This would suggest that improvements in amide hydration free energies in the Drude-2013 model are primarily due to stronger carbonyl–water interactions. Increasing the N—H partial charges and enhancing N polarizabilities were both successful in improving the solubility of the backbone, with the latter also stabilizing the GB1 β-hairpin.

Adjusting parameters in any model should be performed with care. Atomic polarizabilities in the Drude-2013 model, for example, have already been scaled from their gas-phase values by 0.85 to reproduce experimental values for the dielectric constants of pure liquids of pyridine and pyrrole.^21,73 Our proposed increase of 30% for Drude-2013 backbone amide N polarizabilities would increase them by 10.5% over their gas-phase values. While the dielectric constant for pyridine, which does not have an N—H group, was unchanged between gas-phase scaling factors of 1.0 and 0.85, the dielectric constant for pyrrole, which does have an N—H group and should be more representative of backbone amide values, was overestimated by ∼20% with a scaling factor of 1.0,⁷³ suggesting that our polarizability changes would destroy the improvement in bulk properties with the 0.85 scaling factor. However, while increasing the atomic polarizabilities above the 0.85 scaling factor clearly has adverse effects on the dielectric constant of pyrrole, changes in the dipole and quadrupole moments were much less significant, differing by only 2%–4% between 0.85 and 1.0, with the higher scaling factor producing quadrupole moments slightly closer to quantum mechanical calculations.⁷³ Since we only increase the polarizability of protein backbone amide N atoms while keeping all other protein atoms at their 0.85 scaled values, the overestimation of the pyrrole dielectric constant should be limited. In addition, for unfolded states of proteins, the bulk properties of water, which we do not change, should be more important than the bulk properties of the protein backbone. Therefore, overestimation of the pyrrole dielectric constant becomes less important in this context. For the folded states of larger proteins, where the dielectric constant near the center of the protein is much smaller than for bulk water and much more dependent on the electrical properties of the amino acids,⁷⁴ changes in the dielectric constant of the backbone will be more significant than for the small β-hairpin of GB1. Studies on larger proteins will be necessary to determine the effects of polarization on buried residues. Finally, a major advantage of the Drude-2013 model is the avoidance of the overpolarization present in many fixed charge force fields, which often leads to overstabilization of secondary structures. It accomplishes this by introducing dynamic atomic polarization while simultaneously reducing some of the intrinsic polarization. We observed that overpolarization in the Drude-2013 model can also lead to overstabilization of the secondary structure, although this appears to be system dependent as opposite effects were observed for the GB1 β-hairpin and the Ala₁₀ α-helix.

The C36 model greatly overstabilizes the GB1 β-hairpin due to a non-native salt bridge that occurs in the turn region. Elimination of this salt bridge allows C22^* to only slightly overstabilize the β-hairpin. Weaker interactions with the TIP3P water model compared to the SWM4-NDP water model in the Drude-2013 systems along with stronger peptide–peptide hydrogen bonding, as well as a more hydrophobic tryptophan side chain, appear to be the stabilizing factors for C22^*. Ultimately, the C22^* and Drude-2013 models produce roughly equivalent folded fractions that are within or nearly within the experimental range of 30%–50%.^5,6 In the recently released C36m force field, the CMAP potential in the C36 force field was refined to more accurately reproduce protein unfolded states observed in experiments.⁷⁵ However, rescaling our C36 REUS trajectories with the C36m force field parameters according to Eq. (4) produced no discernible difference in the folding free energy landscape of GB1 (see Fig. S12 of the supplementary material), in agreement with the finding that C36 and C36m sample similar conformational ensembles for GB1.⁷⁵ The C36m force field also includes a pair-specific Lennard-Jones parameter to weaken salt bridge interactions between side chain guanidinium N (Arg) and carboxylate O (Asp/Glu) atoms to correct for an overestimation of the equilibrium association constant of a guanidinium-acetate solution and an underestimation of its osmotic pressure.⁷⁵ This was proposed as a simpler, less perturbative method to that employed by Piana et al. in the C22^* force field, which reduced the partial charges of the guanidinium and carboxylate groups to achieve the same effect.²⁴ The additional carboxylate–ammonium-specific Lennard-Jones parameter could improve the folding properties of GB1 in the C36m force field by preventing the formation of the non-native Asp47–Lys50 salt bridge.

Finally, we observe that tryptophan fluorescence, an oft-used method for studying protein folding, does not properly capture the full folding pathway of the GB1 β-hairpin. For all three CHARMM models studied here, hydrophobic collapse precedes hydrogen-bond formation, and the model of hydrophobic contacts proposed by Muñoz et al. does not hold for our simulations.⁹ GB1 folds via a “zippering” mechanism, with hydrogen bonds in the turn region forming first and hydrogen bonds near the termini forming last. In this mechanism, tryptophan fluorescence would underestimate the true folding rate since the energetic barrier for folding is hydrogen bond formation, not a hydrophobic collapse. Another method that monitors non-native contacts is an alternative approach that has been proven to correctly measure folding rates.⁶⁷ Additionally, one could combine tryptophan fluorescence with another technique that also captures the zippering mechanism observed in our simulations, such as FRET to measure contacts near the termini.

See supplementary material for the modification of the partial charges and Lennard-Jones parameters, analysis of the dipole moments, Ala₁₀ folding free energies, and analysis of hydrogen bonding.

This work was supported by a National Science Foundation Award (Grant No. MCB-1452464) to J.C.G. and by NSERC of Canada through the Discovery Grant Program (Application No. 418505-2012) to C.N.R. Computational resources were provided via the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF Grant No. OCI-1053575 and Compute Canada (RAPI: djk-615-ab). E.T.W. thanks NSERC of Canada for an Undergraduate Student Research Award. The authors also thank Benoît Roux for helpful discussions.