Polarizable and non-polarizable force fields: Protein folding, unfolding, and misfolding

The physiological function and physico-chemical properties of biomolecules are inherently linked to their conformational ensembles.^1–3 Computational techniques, ranging from structure refinement tools^4–6 to simulations of protein motions,^7–9 are indispensable for the study of thermodynamics, kinetics, and mechanisms encoded within conformational ensembles.¹⁰ In particular, in the field of biomolecular simulations, they have repeatedly proven to be of great value for understanding and predicting macromolecular quantities based on fundamental interaction potentials.^11,12

The study of how proteins fold into their functional three-dimensional structure, i.e., their native state, is a prominent example in this respect. Based on the idea that the native state is thermodynamically stable at physiological conditions¹³ and on the fact that a random search of the conformational space is incompatible with the experimentally observed folding times,¹⁴ the concept of folding funnels emerged.^15,16 Following this concept, the native state represents the most favorable minimum in a protein’s free energy landscape (FEL). Unfolded conformational states are characterized by higher relative free energies and conformational entropies compared to the native state minimum. Hence, the FEL is said to have the shape of a funnel. It has to be borne in mind that the folding funnel theory represents the system in a simplified way.^17,18 The native state, for example, is not one static structure but characterized by constant conformational rearrangements.¹⁹ However, these motions are small and fast compared to transitions to the unfolded conformational ensemble. Hence, when we refer to the native state, we really discuss the native state ensemble.

Numerical studies are long established in the exploration of protein folding and unfolding processes,^20–24 but the quality of the results hinges on sufficiently accurate descriptions of both the native and unfolded states. This is tied to another central problem in structural biology, i.e., the prediction of a protein’s native state from its sequence.²⁵ The ungraspable complexity of a protein’s FEL renders de novo structure prediction a complex task, but continuous efforts have led to significant methodological advancements.^26,27

While predictions based on sequence homology and rotameric libraries are still the go-to approach, folding predictions based on molecular simulations become ever more relevant. Two important reasons for this are the significant progress in the theoretical (force field based) description of large biomolecular systems and the steady increase in computational power. The majority of current force fields used in biomolecular simulations approximate internal interactions as harmonic or periodic functions, while non-bonded interactions are based on pair-wise additive potentials.^28,29 Since the electrostatic interactions are typically described by point charges located at the atom centers, polarization effects are inherently neglected.

In an effort to improve the description of non-bonded interactions, in particular electrostatics, polarizable force fields have been developed.^30–33 It was shown that the introduction of polarizabilities can alter the relative stability of secondary structure elements^34,35 and their folding and unfolding behavior.^36,37

In this study, we provide a survey of protein folding and unfolding dynamics in three non-polarizable force fields, i.e., AMBER ff14SB,³⁸ CHARMM 36m,³⁹ and GROMOS 54A7,^40,41 and the polarizable Drude 2013 force field.⁴² All four force fields were developed for MD simulations of proteins but differ in their parameterization strategies.^43–45 Most notably, Drude 2013 accounts for the induced electronic polarization inherent to conformational rearrangements through classical Drude oscillators.⁴² While the implementation of Drude 2013 is among the most efficient polarizable force fields in terms of calculation time, it is still approximately four times slower than non-polarizable additive force fields.⁴⁶ We note that very recently, an update in the form of the Drude 2019 force field was presented,⁴⁷ which is not considered in the present study.

AMBER ff14SB has been reported to surpass the preceding AMBER force fields as it improves helical stability and generally provides highly reliable backbone dynamics of folded proteins.³⁸ Similarly, one of the main objectives in the development of GROMOS 54A7 was the correction of a previously observed slight destabilization of α-helical structures.^41,48 CHARMM 36m, on the other hand, was specifically tuned to describe both the folded and disordered conformational ensembles, with a particular focus on intrinsically disordered proteins (IDPs).³⁹ As it became clear that IDPs play a central role in allosteric cellular signaling associated with tumor growth, inflammation, and cardiovascular disease, the need to explore their characteristically shallow free energy landscapes increased.^49–51 

An accurate description of disordered conformational states is also pivotal when studying partial unfolding in well-structured globular proteins, which is often central for their function and biophysical properties.^52–54 However, the detailed structural characterization of disordered conformational ensembles still remains a major challenge⁵⁰ despite substantial experimental and computational efforts in this direction.^55–60 

Another aspect that has to be considered is the role of the water model used as several recent studies report a significant influence of the water model on the behavior of biomolecules.^61–69 Within this study, the aim is to compare different force fields in combination with their respective standard water models, i.e., TIP3P⁷⁰ for AMBER ff14SB, CHARMM-modified TIP3P (TIPS3P)⁷¹ for CHARMM 36m, SWM4-NDP⁷² for Drude 2013, and SPC⁷³ for GROMOS 54A7. Consequently, we can neither compare effects stemming exclusively from the different force fields nor compare effects stemming exclusively from the different water models. For brevity, we only name the different force fields in the following, but we stress that with this, we actually refer to the respective force field–water combination.

To understand the fundamental molecular driving forces governing folding dynamics, typically, fast-folding proteins have been studied as model systems.^79,80 For this survey, we selected three such systems to evaluate differences in the ensembles captured with the four force fields described above. The proteins in question are variants of chignolin (CLN025),⁷⁵ villin,⁷⁷ and a WW domain.⁷⁸ More detailed information on these systems is provided in Table I.

Our main focus lies on the de novo mini-protein CLN025.⁷⁵ Its structure and dynamics have been thoroughly characterized with experiment and long-timescale simulations.^81–83 With only ten residues, CLN025 is the smallest model system in our set, which also shows the fastest folding time (<1 µs).⁷⁴ With this study, we focus mainly on force field based differences in the thermodynamics and kinetics of the CLN025 ensemble. For an overview on the associated folding mechanisms, we refer the reader to a recent work from McKiernan.⁸³ 

In addition, we consider a villin headpiece variant and a WW domain, which represent standard systems for fast-folding α-helical and β-sheet proteins, respectively. The villin headpiece subdomain consists of three α-helices surrounding a compact hydrophobic core (see Fig. 1). As a centerpiece of protein folding studies, villin headpiece variants were studied extensively in experiments^77,84–91 and numerical studies.^{22,74,87,92–96} Folding times obtained in experiments range from 0.7 µs to 4 µs^77,84,85 for villin. The variant considered in the present work is a norleucine/norleucine double mutant containing 35 residues (see Table I),⁷⁷ for which folding times of 0.7 µs (experiment)⁷⁷ and 2.8 µs (simulation)⁷⁴ were reported.

WW domains are relatively small proteins forming anti-parallel β-sheets and contain two eponymous tryptophan residues (see Fig. 1). Multiple variants were studied in experiments exhibiting folding times between 4 µs and 800 µs.^78,97–102 This range of folding times spanning several orders of magnitude was achieved through mutations in the two loops and the replacement of the second tryptophan residue (at the end of the third β-strand) with phenylalanine.^78,98,99,101 Early attempts to study WW domain folding in computer simulations failed,¹⁰³ which is likely related to the used CHARMM22 force field.^104,105 However, later works consistently succeeded in elucidating the folding mechanism.^{74,93,101,106–110} The WW domain variant considered in the present work was only modified in loop L₁ and contains 38 residues (see Table I).⁷⁸ In experiments, the folding time was determined to be 12 µs.⁹⁹ 

For all protein-force field combinations, we performed three independent NPT simulations at 337 K and 1 bar, one simulation starting from the folded PDB (protein databank) structure (F) and two simulations starting from (partially) unfolded states. This choice of temperature is motivated by the early study of Freddolino et al.¹⁰³ To generate the unfolded states, we simulated the solvated systems at 500 K using the CHARMM 36m force field. The thus obtained configurations were subsequently clustered based on C_α-RMSD. For each system, two cluster representatives were picked, which served as unfolded starting structures. These two structures will be named U1 and U2 in the following. U1 has always a lower C_α-RMSD with respect to the crystal structure than U2, i.e., is less unfolded (see Fig. 1). The protonation states of each system are the same for all force fields and different starting structures (see Table I for details).

The parameters for norleucine (two present in the studied villin variant) are only incorporated into CHARMM 36m. For the other three force fields considered here, it had to be parameterized. The AMBER ff14SB parameters were derived using the antechamber module from AmberTools 19¹¹¹ using a capped norleucine in order to avoid charge artifacts. The AM1-BCC charge method was used to derive partial charges in combination with AMBER atom types. Main chain atoms were defined during parameterization in order to avoid sequence dependency artifacts. For Drude 2013 and GROMOS 54A7, the backbone and the side chain parameters up to and including the C_δH₂ group were described using the corresponding parameters from lysine. The C_ϵH₃ group was described using the CH₃ group parameters from isoleucine.

For all simulations with the AMBER package,¹¹¹ topology files and starting coordinates were generated with the program tLEaP of AmberTools 19.¹¹¹ The AMBER ff14SB force field³⁸ was used to describe the proteins, along with the TIP3P water model.⁷⁰ All systems were placed in a cuboidal box with a 10 Å padding distance in each direction. Before production simulations, all systems were subjected to an elaborate equilibration protocol.¹¹² All AMBER simulations were performed with the graphical processing unit (GPU) implementation of the pmemd module of AMBER 18.¹¹¹ A constant temperature of 337 K was maintained with the Langevin thermostat¹¹³ at a collision frequency of 2 ps⁻¹. Constant atmospheric pressure was maintained with the Berendsen barostat¹¹⁴ with a pressure relaxation time of 2 ps. Beyond a non-bonded cutoff of 8 Å, long range electrostatics were treated with a particle mesh Ewald method.¹¹⁵ The charges of the proteins were neutralized by a uniform background charge. All bonds involving hydrogens were restrained using the SHAKE algorithm,¹¹⁶ to allow for a time step of 2 fs. Each simulation was run for a total of 10 µs, collecting frames every 10 ps.

The simulations with the CHARMM 36m force field were conducted with OpenMM 7.3.1.¹¹⁷ The starting structures were placed in a cubic simulation box and solvated in CHARMM-modified TIP3P (TIPS3P) water⁷¹ using a minimal padding distance of 10 Å. The charges of the proteins were neutralized by the addition of ions (see Table I). After energy minimization and a short equilibration (40 ps), the production runs were started (the short simulations were only preformed for the unfolded structures). The equations of motion were integrated using Langevin dynamics with a friction constant of 0.1 ps⁻¹ and a time step of 2 fs. The Coulombic interactions were calculated using the particle mesh Ewald method,¹¹⁸ and a real space cutoff of 12 Å was used. The Lennard-Jones interactions are switched to zero between 10 Å and 12 Å. Pressure coupling was performed with a Monte Carlo barostat.^119,120

The simulations with the Drude 2013 force field were conducted with OpenMM 7.3.1.¹¹⁷ The starting structures were generated from the CHARMM starting structures (following minimization and short equilibration), and TIPS3P water was replaced with polarizable SWM4-NDP water.⁷² The equations of motion were integrated using a dual Langevin integrator based on the extended Lagrangian method^121,122 using a time step of 1 fs. The regular atoms are coupled to a thermostat at 337 K with a friction constant of 0.1 ps⁻¹, while the Drude oscillators are coupled to a thermostat at 1 K with a friction constant of 20 ps⁻¹. A hard wall constraint at 0.2 Å was used. The Lennard-Jones and Coulomb interactions as well as the pressure coupling were handled as in the CHARMM case.

The simulations with the GROMOS 54A7 force field were conducted with GROMACS 2019.2.¹²³ The starting structures were placed in a cubic simulation box and solvated in SPC water⁷³ using a minimal padding distance of 10 Å. The charges of the proteins were neutralized by the addition of ions (see Table I). After energy minimization and short NVT and NPT equilibrations (1 ps each), the production runs were started. The equations of motion were integrated using Langevin dynamics with a friction constant of 0.1 ps⁻¹ and a time step of 2 fs. The Coulombic interactions were calculated using the particle mesh Ewald method.¹¹⁸ An identical cutoff of 12 Å was used for both the Lennard-Jones and real space Coulomb interactions. Pressure coupling was performed with the Parrinello–Rahman barostat,¹²⁴ and bond constraints were maintained using the LINCS algorithm.¹²⁵ 

To analyze the temporal evolution of the systems, we evaluate three quantities for each frame. We assign a secondary structure element to every residue making use of the DSSP algorithm,¹²⁶ determine the fraction of native contacts Q as introduced by Best, Hummer, and Eaton,¹²⁷ and calculate the root-mean-square deviation of all C_α-atoms with respect to the native structures (C_α-RMSD).

We also perform a principal component analysis (PCA) and time-lagged independent component analysis (TICA). For both the analyses, we use inverse distances between the C_α-atoms as input features. For CLN025, we consider all residues in these analyses; for villin and WW domain, however, we only include residues of the structured core, i.e., residues 2–33 for villin and residues 10–34 for the WW domain. We note that we refer to residue 1 as the first residue in the studied crystal structures and refer to the following residues using sequential numbering. For the TICA,^128,129 a lag time of 25 ns was chosen.

For CLN025, we built Markov state models (MSMs) to extract thermodynamic and kinetic information from the trajectories. MSMs are statistical models that describe the dynamical processes in simulations through the transition probabilities between conformational states. As the first step, for the construction of MSMs, the trajectories need to be converted into a succession of discrete states. The transitions between these states are counted with consideration of the lag time, and the transition probabilities are calculated. Here, we clustered the CLN025 trajectories of each force field into 200 clusters using k-means clustering based on the three TICs with the largest eigenvalues. Bayesian MSMs^130,131 were constructed on these discretized trajectories using a lag time of 25 ns. For the Drude 2013 force field, only 63% of the frames are connected and therefore included in the MSM. For each force field, we selected the appropriate number of PCCA+ (Perron-cluster cluster analysis+) clusters¹³² based on the gap between successive eigenvalues to coarse-grain the MSMs. The MSMs of all force fields have been validated with the Chapman–Kolmogorov test.¹³⁰ Stationary probabilities (i.e., equilibrium probabilities) and mean first passage times (i.e., average time until a transition takes place; inverse of the rate) have been calculated from the MSMs. The distributions of native contacts Q and C_α-RMSD values have been weighted with the calculated probabilities from the MSMs and separately for the resulting metastable states according to the PCCA+. To retrieve representative structures for each metastable state, we identify all frames that show a likelihood above 90% to be a member of that state. We then perform a hierarchical clustering on the C_α-positions with a distance cutoff of 0.25 Å. The highest populated structure is then represented as the most prominent in the MSM illustrations.

For villin and the WW domain, we retrieve representative structures from each trajectory based on a hierarchical agglomerative clustering on the C_α-positions applying a distance cutoff of 8 Å for villin and 6 Å for the WW domain. Subsequently, we identify the cluster representative closest to the native structure in terms of C_α-RMSD and depict it as an overlay with the crystallized native conformation.

All described analyses were performed, in part, by using cpptraj and pytraj from AmberTools 19¹¹¹ and the PyEMMA package.¹³³ 

The smallest system considered in this study is the ten-residue β-hairpin CLN025. As described in Sec. III, we study the folding dynamics of CLN025 based on simulations started from the crystallized native state (CLN025-F) and two unfolded conformations (CLN025-U1 and CLN025-U2). Figure 2 shows the time evolution of the secondary structure elements, native contacts Q, and the similarity to the native conformation as C_α-RMSD for the simulations of CLN025-F. These analyses already indicate several folding and unfolding events when using either AMBER ff14SB or CHARMM 36m. While in AMBER ff14SB, the folded state is clearly more prevalent, folded states appear for roughly half the time in the CHARMM 36m trajectory.

Simulations performed with the Drude 2013 force field only appear to explore disordered and misfolded conformational states as even the simulation started from the native structure unfolds within the first few nanoseconds. For t > 9 µs, a rather stable bend segment appears, which is joined by a second bend segment at about 14 µs [see Fig. 2(D)]. Simulations with the GROMOS 54A7 force field fluctuate between the folded and unfolded states in rapid succession. The data obtained from the unfolded starting structures CLN025-U1 and CLN025-U2 yield qualitatively comparable results (see Figs. S1 and S2 of the supplementary material).

We further perform an in-depth analysis of the associated kinetics and thermodynamics using MSMs. Building MSMs for the different force fields allows us to identify differences in state probabilities and transition rates. To this end, we combine the CLN025 simulations that originate from the different starting structures. Additionally, we apply the same workflow a 106 µs trajectory of CLN025,⁷⁴ kindly provided by D.E. Shaw Research (DESRES), as a benchmark for our results. We note that the CLN025 trajectories of the present study were simulated at a slightly lower temperature (337 K) than the DESRES trajectories (340 K) and that the DESRES trajectory was obtained using CHARMM22^*. In Fig. 3, we depict the captured conformational ensemble for each force field as projection onto the DESRES TICA space. Thereby, we can directly compare the resulting folding free energy landscapes. Additionally, we show representative structures of the conformational states identified in the MSM workflow including their state population and mean first passage times (see also Tables S-I and S-II of the supplementary material).

Already from the FEL, we can identify clear distinction between the force fields in our set. AMBER ff14SB and CHARMM 36m show a highly similar pattern to the DESRES reference trajectory, with two distinctly separated minima, which represent the native (folded) and unfolded conformational states. However, the depth of the minima differs significantly between the sampled ensembles. The minima observed in the DESRES reference trajectory translate into a probability of 78% for the folded state. AMBER ff14SB represents the native conformations with a probability of 92%. The CHARMM 36m ensemble, on the other hand, is characterized by a probability of 47% for the folded state and two unfolded states with probabilities of 50% and 3%, respectively. This observed shift of probabilities away from the native state is well in line with CHARMM 36m being specifically parameterized to sample disordered conformational ensembles. We also note that the probabilities of the folded states are consistent with the results shown in Figs. 2(A) and 2(C).

The conformational ensembles and corresponding FELs captured with Drude 2013 and GROMOS 54A7 show substantially different characteristics than the reference DESRES trajectory. However, as described above and highlighted in Fig. 2(D), even when starting from the native conformation, Drude 2013 unfolds the β-hairpin within a few nanoseconds. As indicated above, also the FEL of CLN025 sampled with GROMOS 54A7 implies significant differences to the conformational ensemble of the DESRES reference trajectory. It appears shallower overall, with lower barriers between the most favorable minima, which translates into the corresponding faster transition time scales. With mean first passage times of 100 ns between the folded and highest populated unfolded states, GROMOS 54A7 clearly exhibits the fastest (un)folding dynamics. This trend can already be seen in Fig. 2 as the native contacts form and break rapidly. In terms of state populations, GROMOS 54A7 shows result as CHARMM 36m. The conformational state that resembles the native state most closely is represented with a probability of 45%. The remaining three conformational states identified with the MSM workflow exhibit probabilities of 44%, 9%, and 2%, respectively. We note that the highest populated unfolded state in the GROMOS 54A7 force field is structurally significantly more similar to the native state than for CHARMM 36m (see Figs. S22 and S24 of the supplementary material).

Villin is a widely used model protein to study the dynamics of α-helix formation. As described in Sec. III, we again perform simulations starting from the native state (VIL-F) and two unfolded conformations (VIL-U1 and VIL-U2) for each force field. The simulations starting from VIL-F are stable in the AMBER ff14SB force field over the studied 10 µs, while in the CHARMM 36m trajectory, which is more than twice as long, the folded structure partially unfolds at t > 12 µs. When villin is studied in the Drude 2013 force field, the folded state is partially retained for the 11 µs studied. In the GROMOS 54A7 force field, villin transitions to a misfolded state.

The unfolding events observed in CHARMM 36m and GROMOS 54A7 indicate a destabilization of the native state since the experimental melting temperature of the WW domain is 23 K higher than the temperature studied here (see Table I). Nevertheless, it appears that in both force fields, helix α₃ (see Fig. 1) is the most stable helix, a finding which is consistent with isotope-edited IR spectroscopy data⁹¹ (a more detailed discussion of the VIL-F simulations can be found in the supplementary material).

The results for the simulations started from the unfolded villin conformation VIL-U1 are shown in Fig. 4. The simulation performed using AMBER ff14SB refolds within 2 µs and then remains close to the native state, with the exception of a partial unfolding event after around 9 µs. The conformational ensemble resulting from simulations starting from VIL-U1 using CHARMM 36m does not include the native state. After numerous conformational rearrangements in the first 15 µs, the system appears to become trapped in a misfolded state for the following 8 µs. In this misfolded state, helices α₂ and α₃ are correctly formed, but instead of helix α₁, a β-sheet is formed. Starting a simulation from the same structure with the Drude 2013 force field shows partial refolding within the first few nanoseconds. However, the simulation does not sample the native state. Neither do the simulations with the GROMOS 54A7 force field. Nevertheless, the analyses shown in Fig. 4(G) indicate numerous partial folding and unfolding events of all three helices. It appears that in simulations using the GROMOS 54A7 force field, helix α₃ is the most likely helix to be folded.

The unfolded villin structure VIL-U2 has a higher C_α-RMSD to the crystal structure than VIL-U1. In Fig. 5, we profile the degree of folding captured with each force field on starting from VIL-U2. Again, the simulation performed using the AMBER ff14SB force field resembles the characteristics of the native state closely in the captured 10 µs. Only slight deviations in helix α₃ remain [see Fig. 5(A)]. For CHARMM 36m, the results are quite different as in this case, the simulation explores helical and β-sheet conformations, but it does not reach the native state. The simulations started from VIL-U2 using Drude 2013 misfolds toward one large α-helix within the initial 2 µs, which spans the residues of helices α₁ and α₂ in the native structure. As the simulation progresses, this helix splits into two, and also helix α₃ is partially formed. Despite not reaching the native state, the system reaches a conformation that is comparable to structures captured with Drude 2013 simulations started from the VIL-F structure (see Fig. S26 of the supplementary material). As shown in Fig. 5(G), simulations with GROMOS 54A7 started from VIL-U2 approach the characteristics of the native conformations. While helices α₁ and α₃ are stably folded throughout most of the 10 µs trajectory, the second helix is only formed in the last microsecond of the trajectory. These last conformations then exhibit characteristics resembling the native state.

Apart from the simulations with CHARMM 36m, all simulations of VIL-U2 at least approach the native state. However, the order in which the helices are formed is different. In AMBER ff14SB, we observe α₁ → α₂ → α₃; in Drude 2013, we observe α₃ → α₁ → α₂ (note that α₁ is not folded correctly); in GROMOS 54A7, we observe α₁ → α₃ → α₂ (note that α₁ unfolds and refolds several times). Experimental results suggest that predominant folding pathways start with the formation of helix α₃,⁹¹ which is observed here only for Drude 2013, but Drude 2013 does not reach the fully folded state within the ≈5 µs trajectory studied. In both AMBER ff14SB and GROMOS 54A7, helix α₁ is formed fist. This discrepancy and the agreement of Drude 2013 with the experimental findings should, however, not be overrated since we here study folding of a single unfolded conformation and not of an ensemble of unfolded conformations as done in experiments.

The statistical significance of the achieved simulation data unfortunately does not allow a reliable analysis of the captured ensembles with MSMs. However, in Fig. 6, we depict the captured conformational ensembles as projections onto the PCA space of a 125 µs DESRES reference trajectory (PC1 explaining 18% of the variance, PC2 11%). This representation allows a direct comparison of the FEL resulting from each starting structure and force field.

All AMBER ff14SB simulations represent the native state as the most favorable minimum, regardless of the starting structure. The most extended starting structure (VIL-U2) also captures a significant region of the unfolded portion of villin’s conformational ensemble. Interestingly, the simulation starting from the partially unfolded conformation (VIL-U1) hardly samples disordered areas of the FEL.

With CHARMM 36m, the unfolded conformational space is represented more prominently than with AMBER ff14SB. Even on starting from the native state, CHARMM 36m captures an extensive portion of the unfolded region. While the simulation starting from the unfolded conformation VIL-U2 explores extensive regions of villin’s FEL, we find no probability density in the proximity of the minimum representing the native state. Starting from the unfolded conformation VIL-U1, on the other hand, results in a FEL that closely resembles the reference DESRES trajectory and also includes the native conformational state.

The FELs of villin as described with Drude 2013 are highly consistent and appear to be rather independent from the starting conformation. However, the minimum sampled in all Drude 2013 simulations does not correspond to the native state. Simulations performed using GROMOS 54A7 lead to very similar observations. The minimum representing the native state is hardly sampled.

To visualize the structural characteristics of each simulation, we clustered each trajectory based on the C_α positions (see Sec. IV). In Fig. 7, we depict the representative structure with the smallest C_α-RMSD to the native conformation for each villin simulation. These representative conformations were also further projected onto the PCA space in Fig. 6.

The studied WW domain consists of 38 residues arranged in three β-strands. While similar in length to villin, we consider its computational study more challenging as the formation of β-sheets is significantly slower than for α-helices.¹³⁴ We characterize the folding dynamics of the WW domain with each force field based on simulations starting from the native state (WW-F), a partially unfolded conformation (WW-U1), and a fully unfolded conformation (WW-U2).

The native WW-F structure is stable in both the AMBER ff14SB and CHARMM 36m simulations over the entire trajectory (10 µs and 20 µs, respectively). In contrast, the simulation performed with Drude 2013 leads to partial unfolding of the protein. Only the central region of the WW domain stays intact. The WW-F structure also appears to be unstable when simulated using the GROMOS 54A7 force field. As in the case of villin, unfolding is not expected as the melting temperature of the WW domain is 13 K higher than the temperature studied here (see Table I). However, in contrast to villin, the WW domain is stable for 20 µs in CHARMM 36m (a more detailed discussion of the WW-F simulations can be found in the supplementary material).

The results of the simulations starting from the partially unfolded WW-U1 structure are depicted in Fig. 8. Both the AMBER ff14SB and CHARMM 36m force fields direct the system back to the native state within the first 2 µs, which is stable for the remaining trajectory. When performing simulations starting from WW-U1 using Drude 2013, we observe further unfolding. Already in the first stages of the simulation, the two main β-strands transition toward disordered conformational states. At ≈3 µs, parts of this region refold as indicated by all three quantities shown in Fig. 8(D). However, at the same time, the residues of the short β-strand (β₃) form an α-helix. At ≈9 µs, this α-helical structure replaces the two main β-strands, resulting in a misfolded state. The dynamics of the GROMOS 54A7 simulations are less severe. While the partially unfolded WW-U1 conformation does not fold to the native state within the ≈10 µs studied, it also does not completely unfold. That is, two of the three β-strands of the WW domain stay partially folded over the whole trajectory. We also observe that both the native contacts and the C_α-RMSD show little fluctuation for t ≥ 5 µs. This indicates that this trajectory is exploring one minimum in the free energy landscape from which it does not escape within the captured timescale.

None of the simulations started from the fully unfolded WW-U2 structure capture the native state of the system within 6 µs–20 µs of simulation time. In the simulation with AMBER ff14SB, the protein misfolds into a stable conformation consisting of short 3₁₀-helices and turns within the first microsecond. For CHARMM 36m, the analyzed quantities indicate that the system does not exhibit a stable fold within the simulated 20 µs. In the Drude 2013 simulation, we observe a partial formation of the two main β-strands. However, this onset of folding does not lead to sampling of the native state conformation. In addition, GROMOS 54A7 does not lead to the native fold but, instead, captures a stable misfolded state including an α-helix (a more detailed discussion of the WW-U2 simulations can be found in the supplementary material).

As these results already imply, we were not able to construct statistically stable MSMs for the WW domain. Hence, in Fig. 9, we represent the FEL explored with each starting conformation and force field individually as a projection onto a combined PCA space, where PC1 comprises 45% of the variance and PC2 13%. (We do not project onto the PCA space given by the DESRES trajectory recorded for the WW domain because a significantly different variant, i.e., GTT, was studied in Ref. 74.) This analysis emphasizes the substantial impact of the WW domain starting structure on the resulting ensemble. As already suggested in Fig. 8, the FEL described with AMBER ff14SB shows that simulations are strongly stabilized once they are close to native state. However, the simulation starting from the WW-U2 conformation remains in a narrow distinct minimum far away from the folded conformational state. CHARMM 36m also consistently samples the narrow region around the minimum representing the native state not only for the simulation starting there (WW-F) but also for the simulation starting from the partially unfolded conformation (WW-U1). Most intriguingly, CHARMM 36m explores remarkably extensive portions of the FEL when started from the fully unfolded state (WW-U2) including a minimum in close vicinity to the native state.

The FEL resulting from Drude 2013 simulations of the native state identifies a minimum in the vicinity of the native state as the most favorable, while the native state itself does not correspond to a minimum. When starting from the partially unfolded conformation, Drude 2013 does not visit the minimum close to the native state but rather explores unfolded conformations. Similar FEL regions are explored in the WW-U2 simulations. GROMOS 54A7 based simulations display a very different behavior based on the constructed FEL in PCA space. The simulation starting from the native state correctly represents the native minimum but also leads to considerable exploration of unfolded conformations. The native minimum is also sampled if the simulations are started from the partially unfolded conformation. However, the simulation starting from the unfolded conformation (WW-U2) does not visit the vicinity of the native state minimum.

The structural characteristics discussed for each simulation are further illustrated in Fig. 10. We again clustered each simulation and depict the representative conformation, which resembles the native conformation most closely (see Sec. IV). The respective conformations are also depicted as a projection onto the PCA space shown in Fig. 9.

In this study, we characterize the folding dynamics of the three fast-folding proteins using four different force fields, i.e., AMBER ff14SB, CHARMM 36m, GROMOS 54A7, and the polarizable force field Drude 2013 in conjunction with their standard water models (TIP3P, TIPS3P, SPC, and SWM4-NDP, respectively). Focusing on the smallest system in our test set, CLN025, we observe several clear differences in the conformational ensembles captured by the different force fields. Due to its comparably small size and fast transitions between the native and disordered conformational states, CLN025 is a prime model system for quantitative studies of β-hairpin folding dynamics. Here, we evaluate the performance of each force field based on simulations starting from the native and two unfolded conformations. Analyses based on the secondary structure, native contacts, and C_α-RMSD of each simulation already suggest that the captured dynamics are independent of the starting structure. While the ensembles do not appear to be biased by the initial conformation, we find significant differences between the considered force fields. To ensure statistically and physically meaningful analyses, we characterize the ensembles resulting from each force field using MSMs. We apply the same MSM workflow to an already published 106 µs DESRES trajectory of CLN025 (simulated with CHARMM22^*) to provide a reference for our evaluations. The MSM built from this ultralong DESRES trajectory assigns a probability of 78% to the folded state at 340 K, and the corresponding work of Lindorff-Larsen et al.⁷⁴ suggests a melting temperature of 381 K for CLN025. Experimental studies, however, estimate a population of ∼55% at 337 K for the native state⁷⁵ and determine a melting temperature of ∼344 K⁷⁶ (see Table I). An overestimation of the folded fraction of CLN025 has also been observed for other force fields.^120,135

Similarly, also the CLN025 ensemble captured with AMBER ff14SB represents the folded state as clearly prevalent. Furthermore, the observed kinetics agree remarkably well with the reference DESRES MSM. A vague experimental reference for the mean first passage times can be derived from T-jump experiments, which report relaxation times of $∼$100 ns⁷⁶ at 328 K. Although the experimental relaxation times were determined 10 K lower than the present simulations, they still indicate that the folding mean first passage times suggested by AMBER ff14SB are slightly more realistic than the ones suggested by the DESRES reference.

CHARMM 36m explores a similar FEL as the AMBER ff14SB simulation, but it estimates a significantly lower population for the native state (47%). Furthermore, the transition rate toward the native state is estimated about an order of magnitude slower in this force field. Hence, the CHARMM 36m MSM retrieves state populations closer to the experimental estimates, while the kinetics are further off than for AMBER ff14SB and DESRES simulations. Furthermore, we find that the CHARMM 36m ensemble does not agree with the DESRES reference as AMBER ff14SB as well. This might appear counterintuitive as the DESRES simulations are based on a force field from the CHARMM family. However, the aim in the development of the studied CHARMM 36m was to provide a force field that would be applicable for both structured and disordered proteins. It is thus not surprising but rather encouraging that this force field is biased toward the exploration of unfolded conformations. Similar to the CHARMM 36m simulations, also the GROMOS 54A7 ensembles agree well with experimental thermodynamics, assigning a 45% probability to the native state. The CLN025 folding dynamics captured with GROMOS 54A7, however, also excel in terms of the estimated kinetics. Here, the transition rates between the native and unfolded conformational states are substantially faster. While DESRES, AMBER ff14SB, and CHARMM 36m estimate mean first passage times between 1 µs and 2 µs, GROMOS 54A7 fluctuates between the native and unfolded conformations within 100 ns. Drude 2013 can also be considered as a member of the CHARMM family and stands out as the only polarizable force field in our test set. Nevertheless, Drude 2013 is also the only force field that does not capture conformational rearrangements toward the native state. Although one of the performed Drude 2013 simulations is started from the native conformation, the structure unfolds rapidly. As no refolding events are observed, the initial frames depicting the native state are not considered in the MSM.

In addition, we study helix folding dynamics based on simulations of villin using polarizable and non-polarizable force fields. Besides the difference in the secondary structure, villin is more than 20 residues larger than CLN025 and thus poses a more complex sampling challenge to the selected force fields. Again, we perform simulations starting from the native state and two unfolded conformations to evaluate simulation stability and folding dynamics. The individual analysis of each simulation highlights AMBER ff14SB as a performance champion in terms of folding toward the native conformation and stable sampling thereof. While all force fields sample the minimum around the native structure when started from there, AMBER ff14SB is the only one that does not exhibit notable indications of unfolding. We, nevertheless, emphasize that the CHARMM 36m simulations are more than twice as long (>20 µs) as the AMBER ff14SB simulations (10 µs), and also in CHARMM 36m, villin is entirely stable in its native state up to t ≈ 12 µs. When we evaluate the simulations started from the unfolded conformations, it is again AMBER ff14SB that refolds the structure closest toward the native state minimum. This result is somewhat surprising, as also GROMOS 54A7 was specifically optimized to stabilize helical structures. Additionally, we had expected Drude 2013 to be particularly well suited for folding villin as polarization plays a key role in helix formation.¹³⁶ 

The most demanding system in our set is a WW domain, which we selected to study the formation of β-sheets. Although this WW domain is comparable to villin in size, it is well established that folding of β-sheets is significantly slower than helix formation.¹³⁴ We analyze simulations starting from the native, a partially and a fully unfolded structure. AMBER ff14SB and CHARMM 36m show highly stable simulations of the native state. Furthermore, both force fields also rearrange the partially unfolded starting structure toward the native conformation within the first 2 µs. Both GROMOS 54A7 and Drude 2013, on the other hand, display lower stability, i.e., partial to full unfolding, in the simulation starting from the native state. For the simulations started from the partially unfolded conformation, the system appears to be trapped in a misfolded conformational state in both force fields. For Drude 2013, this misfolded state even includes an α-helix. When starting from the fully unfolded conformation, none of the considered force fields is able to refold the WW domain to its native structure. In all runs except Drude 2013, α-helical structures appear. While it has been reported before that β-sheet formation can proceed via α-helical intermediates,^137,138 only the simulation with CHARMM 36m shows a direct transition from an α-helical arrangement to an almost folded WW domain, which, however, is short lived (see Fig. S32 of the supplementary material). The longer lived α-helical structures found in AMBER ff14SB and GROMOS 54A7 point to an over-stabilization of this structural feature.

In summary, based on the simulation setup, model systems, and analyses presented in this study, we report the following key observations. CHARMM 36m performs slightly better than AMBER ff14SB in sampling representative state probabilities, while AMBER ff14SB performs slightly better in α-helix formation. However, we note that AMBER ff14SB also quickly misfolds the fully unstructured conformation of the studied WW domain into a highly stable helical conformation. This could be an indication that the helical structure elements are overstabilized within this force field. It has been described for preceding force fields in the AMBER family that the helical propensities of amino acids are generally overestimated.^139,140

The preference for α-helical structures in AMBER ff14SB observed here is likely a feature related to the standard water model of AMBER ff14SB TIP3P. Past studies found that TIP3P induces overly compact protein structures^65,67,141 and that AMBER ff14SB per se underestimates helical propensities.⁶⁹ Generally, it has been shown that a change in the water model can have a significant impact on folding equilibria and kinetics.^61,66,67,69 For instance, it was recently suggested that proteins can refold at very low temperatures,¹⁴² a phenomenon that is likely related to the complex behavior of supercooled water^143,144 and hence strongly depends on the quality of the employed water model. In this survey, however, we always used the respective standard water models of the four force fields considered. Hence, we are only able to compare these four different force field–water combinations and cannot draw conclusions regarding the influence of water models on the protein folding dynamics observed.

While we expected the most accurate results for Drude 2013 simulations, we consistently find it to be outperformed by the other three force fields. Considering its sophisticated treatment of electrostatics and the resulting significant computational overhead, we present simulations of notable length. Nevertheless, we observe surprisingly unstable sampling for all the systems starting from the native conformation. Indications that Drude 2013 destabilizes β-hairpins were already noted in the past,³⁵ and during the preparation of this article, the Drude developers published a work with similar findings.⁴⁷ In this work, Lin et al.⁴⁷ also present an optimized version of the Drude force field (Drude 2019), which is expected to resolve the observed issues. In fact, it was shown that the folded state of CLN025 is stable in Drude 2019 throughout a 1 µs trajectory simulated at 300 K.⁴⁷ 

In the presented work, we profile folding dynamics captured with the polarizable Drude 2013 force field and three non-polarizable force fields, i.e., AMBER ff14SB, CHARMM 36m, and GROMOS 54A7 in conjunction with their respective standard water models. In performing long-timescale simulations of the three model systems, we evaluate the performance of each force field in sampling β-hairpin, α-helical, and β-sheet structures. In particular, for Drude 2013, protein dynamics on such extended timescales have not been studied before. Despite the extensive sampling and sophisticated description of electrostatics, we consistently observe destabilization of the native state ensemble in all Drude 2013 simulations, independently of its respective secondary structure elements. However, the recently published Drude 2019 force field was explicitly developed to address and resolve this issue.⁴⁷ All non-polarizable force fields in our test set perform well in sampling the dynamics within the folded state ensemble. For the ultra-fast-folding mini-protein CLN025, ensembles from both the DESRES and the AMBER ff14SB simulations overestimate the population of the native state. This common limitation of protein force fields is overcome by GROMOS 54A7 and CHARMM 36m, which was particularly developed to sample structured and disordered proteins. Our results further indicate that GROMOS 54A7 even outperforms CHARMM 36m in terms of accurate folding kinetics of CLN025.

In the description of helix formation, AMBER ff14SB achieves the most promising results as even the most extended starting conformation of villin closely folds toward the native state within 10 µs. However, considering our findings for the WW domain, it appears that helices are generally overstabilized in our AMBER ff14SB simulations.

Overall, our findings suggest that current additive force fields are well suited to study not only ultra-fast-folding dynamics but also partial unfolding and refolding processes in feasible computational time.

See the supplementary material for additional figures, tables, and a corresponding discussion.

A.S.K. and P.H.H. contributed equally to this work.

All authors thank the Austrian Science Fund (FWF) for financial support (Project Nos. P30565 and P30737). P.H.H. also acknowledges funding by the University of Innsbruck (NWF-Project No. 282396). The computational results presented were achieved, in part, using the HPC infrastructure LEO of the University of Innsbruck and the Vienna Scientific Cluster (VSC). They are grateful to D. E. Shaw Research for supplying their trajectories of the systems considered in this study.

The data that support the findings of this study are available from the corresponding author upon reasonable request.