Evolution of all-atom protein force fields to improve local and global properties

Introduction

Studies of protein structure and function have historically focused on folded proteins or protein domains, in a large part because they were amenable to experimental structure determination. However, it is now well appreciated that proteins which are either partially or fully disordered perform vital functions .^1,2 In contrast to folded proteins, these so-called intrinsically disordered proteins (IDPs), lacking a well-defined three-dimensional structure, can adopt a broad distribution of conformations^3,4 and cannot easily be studied by conventional structural biology methods. Only a limited number of experimental methods are applicable, most notably nuclear magnetic resonance (NMR),⁵ Förster resonance energy transfer (FRET)⁶ and small-angle X-ray scattering (SAXS).⁷ In addition, experimental IDP data sets are typically averaged over a broad conformational ensemble,^3,4 increasing the difficulty of extracting the underlying structural properties.^8–12

In this context, molecular simulation is clearly a useful tool, as it can potentially yield information on the structural and dynamic properties of IDPs associated with their functions with high spatiotemporal resolution .^13–15 Experimental data may also be incorporated in molecular simulations, either as a bias or to reweight the configurations sampled in an existing trajectory.^8,16–20 Regardless of whether or not such a bias is used, having an accurate simulation force field is expected to improve the quality of the resulting structural ensemble.^{14,15,21–23}

However, as has been well documented in recent years, many force fields in current usage have significant shortcomings when applied to unfolded and disordered proteins ,^22,24–26 often being too compact or giving a poor reproduction of partial secondary structure content in unfolded states (reviewed in^27–30). This is not because they were explicitly optimized to reproduce folded protein structures: in fact, the parameters were usually optimized to describe the properties of small molecule fragments containing functional groups common in proteins, and then assumed to be mostly “transferable” to the macromolecular force field^31–34 although more recent efforts at backbone parameterization have included data derived from the Protein Data Bank.^35–37 Nonetheless, the focus of the systems used in testing and in applications was on folded states. As a result of some 40 years of optimization, all-atom force fields have now been shown to be of sufficient quality to reversibly fold many small proteins (less than 100 residues), a strong validation of the optimization strategy used,^38–41notwithstanding the remaining deficiencies for unfolded states, mentioned above.

Disordered proteins populate a diverse ensemble of conformers of similar free energy and their population is sensitive to small errors in the energy function used for all-atom simulations. Therefore, ensemble-averaged experimental data for IDPs and weakly structured peptides in different solvent conditions are a good target to optimize force fields. Two major shortcomings of force fields with respect to reproducing IDP properties have been identified: incorrect prediction of secondary structure and a tendency to be too collapsed with respect to experiment. The former problem had long been recognized, leading in the case of the Amber force fields to a succession of changes in which the secondary structure propensity was altered by modifying torsion angle parameters (Amber ff94 to Amber ff96, Amber ff99, Amber ff99SB and beyond).^42–48 With the development of advanced sampling techniques and improved computer power, it became practical to use equilibrium data for weakly structured peptides in order to optimize secondary structure propensity, as was done in the Amber ff03*,⁴⁹ Amber ff99SB*,⁴⁹ CHARMM 22*.⁵⁰ and CHARMM 36⁵¹ force fields, to give just a few examples. As a consequence of these and related efforts, force fields can now more accurately capture the propensity for local, α, β or extended structures in disordered proteins.

More recently, it was found that the radius of gyration and related properties of disordered proteins were poorly reproduced by current force-fields: simulations of unfolded or disordered proteins frequently resulted in collapsed structures, often barely larger than folded states,^24,52–54 pointing to insufficiently favorable solvation, a picture that was supported by overly favorable protein-protein interactions in several simulation force fields.^55,56Initial attempts to correct this by changing from the widely-used TIP3P water model⁵⁷ to more accurate models such as TIP4P-Ew^25,52 or TIP4P/2005⁵⁸ yielded some improvements, but disordered states were still too collapsed.⁵²

Two main approaches have been taken to address the issue of collapsed unfolded states. In the first, the strength of the protein-water Lennard-Jones term^22,25 is directly modified, while in the second the water model itself is reparametrized to have a stronger contribution from dispersion interactions.⁵⁹ In both of these instances, and other related studies,⁶⁰ the major change from previous force fields has been to make protein solvation more favorable. These efforts have resulted in significant improvements in capturing not only dimensions of unfolded protein but also intrachain peptide dynamics in disordered peptides.⁶¹ These models have been used extensively to obtain molecular insights into the equilibrium structural ensemble of disordered proteins which otherwise is challenging to obtain solely from experimental data.^13,14,27 In addition to the above efforts, a host of recent studies have been directed at further improvements in reproducing the properties of both folded and disordered proteins in all-atom simulations.^{15,23,62–66}

Because of the known deficiencies of force fields, many investigators now compare their simulation results to experimental data in order to identify potential problems. However, as we will show, agreement with one set of data is not necessarily a strong validation of the results, as different observables can reflect different properties of the ensemble. In order to systematically address this issue, we focus on comparing a series of sequential modifications of the Amber ff03 force field⁴⁵ which have been developed in our laboratories, and for which we have a large set of data available. Using a large data set of disordered and unfolded proteins, we provide a comprehensive survey of Amber ff03*, ff03w, and ff03ws, the evolution of which will be detailed in the next section. We find that the degree of compaction is most significantly affected over this series, with the chain dimensions from Amber ff03ws being most similar to experimental estimates.

Force Field Evolution

The modified versions of the Amber ff03 force field considered here are ff03*, ff03w, and ff03ws. The parent force field ff03 is an all-atom molecular mechanical potential based on the original Amber ff94 force field,³¹ but with atomic partial charges refitted to an electrostatic potential calculated from a quantum chemistry calculation including an implicit solvation model, as well as refitted dihedral potentials.⁴⁵ The model is coupled with a 3-site water model, TIP3P,⁵⁷ to represent the proteins in an explicit aqueous solution. Over the last decade, our groups have developed a number of cumulative improvements to this force field. Prior to the introduction of empirical backbone torsional correction term,⁴⁹ secondary structure bias was a well-known problem. The backbone correction term adjusted for the very small, but nonetheless significant, biases towards one or other secondary structure. A single additional torsion term was applied to the Ψ torsion angle and its parameters optimized to match helix populations inferred from chemical shift data for the (AAQAA)₃ peptide.⁴⁹ This simple correction has been made in the context of the original force fields Amber ff99SB, ff03 and CHARMM22, yielding respectively the ff99SB* and ff03*⁴⁹ and CHARMM22*³⁹ and have been shown to allow the use of same force field to fold both helix and β-strand-dominant peptides.^39–41

Protein-water interactions and water-mediated interactions are clearly important driving forces behind biological assembly processes. Therefore, both protein-water and water-water interactions are critical components in the modeling of biomolecules in aqueous solutions. Explicit modeling of water has been area of intense research activity which has led to much improved models .^67–69 We have focused on the TIP4P/2005 water model,⁶⁹ one of the most accurate four-site point-charge models for liquid water. As a first step, this water model was combined with Amber ff03, together with a backbone correction for secondary structure propensity, resulting in the ff03w force field (a similar approach was applied to Amber ff99SB to give ff99SBw).⁵⁸ This combination yielded some immediate improvements in the thermodynamics of protein folding as well as aqueous solution properties of IDPs such as the temperature-dependent collapse seen in SAXS and FRET experiments.^70–72

Despite the improvements obtained through the backbone torsion corrections and the use of better water models such as TIP4P/2005 and TI4P-Ew, it was found that unfolded and disordered proteins were still overly compact compared to experimental estimates ,^{24–26,52,53,58} suggesting that (in simulations) water acts as a poor solvent for polypeptides.

In an effort to remedy the poor solvation of proteins, we proposed a simple and minimal force field change, by introducing a general scaling factor included in the geometric combination rule for the Lennard-Jones ε parameters between the oxygen atoms in water and every types of atoms in protein.²² Note that in principle an amino acid-based, or atom-based scaling factor could be introduced,²⁵ but would require more experimental data for optimization. We found that for Amber ff03w, a scaling factor of 1.10 (10% increase in the protein-water interactions) has been found to be optimal for matching the simulated and experimental R_g of a fragment of CspTm. Amber ff03w with this scaling of protein-water interactions was named Amber ff03ws.²² The same change in the context of Amber ff99SB-ildn-Q⁷³ was termed Amber ff99SBws. Amber ff03ws has then been shown to reproduce experimentally measured sizes of a number of unfolded or disordered proteins and resolve the nonspecific protein-protein interactions problem.^14,22,26 For compatibility with the new Amber ff03ws, we have also optimized protein-denaturant interactions based on the KBFF denaturant force field^74,75 and extended Amber ff03ws to work with both urea and GdmCl (i.e. KBFFs).⁷⁶

In this work, we have used the Amber ff03-series of force field modifications described above to sample a broad variety of proteins, spanning a range of sizes and sequence characteristics. We will discuss how each of these successive improvements to Amber ff03 family affect the local and global structural properties of IDPs in the next two sections. Methodological details, such as the advanced sampling methods, algorithms, ensembles, and analysis are provided in Supporting Information. A summary of the properties of all the sequences in this work are presented in Figure S1 and simulation run lengths are summarized in Table S1–S4.

In general, any polymer is able to undergo a collapse transition as conditions such as temperature or solvent are varied, from a state in which it forms a completely collapsed globule and attractive interactions dominate, to one in which repulsive interactions dominate and the chain is expanded, the self-avoiding walk limit. At a point in between these extremes, known as the θ state, the attractive and repulsive interactions effectively cancel, yielding chain statistics similar to those of a ghost chain.⁷⁸ Thus, the degree of collapse and its response to changes in environmental conditions such as temperature or solvents is an important global property of intrinsically disordered proteins.

The degree of collapse of a polymer chain can be characterized using a scaling law of the form⁷⁹

in which R_g is the radius of gyration, k is a prefecture, ν is a scaling exponent, and N is the chain length (here we follow the convention in the protein biophysics literature⁸⁰ that N is the number of residues rather than the number of bonds N −1; in practice, this makes a negligible difference). By measuring the internal distance scaling within different proteins, Schuler and coworkers determined average scaling exponents of 0.46 and 0.62 in aqueous solution and in high denaturant solution, respectively, by using FRET experiments,⁶ although there is more variation in water than in denaturant .^6,81 Similarly, Plaxo and coworkers have reported a scaling exponent of about 0.6 when fitting R_g from SAXS experiments under denaturing conditions.⁸⁰ More recent analysis of SAXS and FRET data confirmed these conclusions, although the estimated scaling exponents in water were slightly larger.^9,12,14,82 Thus, proteins in water are close to the θ-state, which is characterized by scaling exponent ν = 0.5 while in denaturant they are closer to a self-avoiding walk, for which ν= 0.6.

To address the scaling behavior in simulations, we compiled simulation data for proteins with a broad range of characteristics. In Figure S1, we list the hydrophobicity, net charge, and praline content for each sequence. The proteins are either intrinsically disordered, or for foldable proteins, in which case we analyze only their equilibrium unfolded ensembles. For small peptides where a folding equilibrium was sampled in the simulations (such as GB1 and TrpCage), we decide the unfolded population according to a cutup using a standard metric such as C_α RMSD or backbone d_RMS. For other foldable proteins (such as CSP, NTL9, Ubiquitin, Barstar, etc.) we did not observe folding transitions within the time scales of the simulation. In Figure 1A we compile the average radii of gyration as a function of the chain length for each protein, and with each force field variant, at 300 K. The average R_g’s approximately follow a power law dependence (Eq. 1) with the fitted parameters and statics tabulated in Table 1. The older generations of the force field, ff03* and ff03w, yield fairly compact ensembles which are much more compact than the experimentally observed unfolded or disordered protein size with scaling exponents of 0.38 and 0.39 respectively, close to those of a compact globule, inconsistent with experiment which shows the chains to be more expanded in water. By contrast, the scaling exponent in water for the latest generation, ff03ws, is much closer to the experimental estimates, with ν = 0.53. In addition, the scaling exponent for ff03ws in denaturant is close to the theoretical value for a self-avoiding walk (due to the finite length of the chain it is slightly larger than that infinite-chain limit). While the trend approximately follows a power law, there are small deviations from the perfect fit, presumably reflecting the broad range of sequence characteristics represented by these peptides (Figure 1A). Notably, the scatter is larger for proteins in water than in denaturant, consistent with a picture in which the denaturant effectively cancels the sequence-specific attractive features of each chain; a similar picture has been observed in experiments .^6,81 We note that we exclude several proteins from Figure 1A as extreme outliers due to their unusually high charge or proline content (Figure S1). These include PROT-C (03w), PROT-N (03w) and tau_174−183 (03*). Proline has computationally been shown to have an impact on the dimension of the disordered proteins in a recent study.⁸³ An all-inclusive version of Fig. 1A is presented in Figure S2.

We have also calculated internal distance scaling separately for each protein, by computing the root mean square distance R_i,j between the C β carbon of every pair of residues(i, j) (where i ≠ j), as a function of sequence separation |i - j|. The scaling of the averaged distances as a function of |i − j| can be approximated by a power law

with b and ∨ the fitting parameters, providing an independent estimate of compactness of the protein.^78,84 Individual fits of the power law function to the scaling of intra-chain distances with sequence separation are presented in Figures S3, S4 and S5, and the distribution of the scaling exponents obtained for each force field is in Figure 1B. The exponent of 0.5 corresponds to a so-called θ -state of the chain in which the effective interactions between monomers are on average neither repulsive nor attractive. Exponents smaller than 0.5 indicate a net attractive interaction between monomers, therefore, a more compact overall conformation. When the interactions between monomers become repulsive, the scaling exponent approaches the limiting value of Ȉ 0.6 for long chains .⁷⁸ Consistent with the scaling of the R_g of the proteins, the average of these exponents is significantly increased when protein-water interactions were optimized in ff03ws. Similarly, there is only a slight increase in ν moving from ff03* to ff03w. This confirms that using a water model which more accurately reproduces the properties of water is not sufficient to improve the accuracy of the all-atom protein simulations in explicit water.

The improvement in the reproduction of the sizes of unfolded and disordered proteins by the newest generation of the force field (ff03ws) is striking. We emphasize that this improvement is not limited to the ensemble-averaged global structural quantities. In earlier work, we have investigated the diffusionand reaction-controlled rates of contact formation within disordered peptides, finding that simulations with the newest generation of the force field more accurately reproduce the experimentally determined dynamics.^61,76 On the other hand, the rates obtained from the older generations of the force field deviate from the experimental data by up to an order of magnitude. We have rationalized this improvement as due to the

We selected three of larger unfolded/disordered proteins in our data set, with a partial secondary structure content, to investigate trends in secondary structure across the series ff03*-ff03w-ff03ws. While the Islet Amyloid PolyPeptide (IAPP)⁸⁵ and the 310–350 fragment from a TAR-DNA binding protein (TDP-43_310−350)¹³ are intrinsically disordered, cold-shock protein (CSP) from Thermotoga maritama⁷² can fold into a five-stranded β -barrel, but only its unfolded state has been sampled in our simulations. We calculate the secondary structure content of the proteins following the widely used DSSP algorithm.⁸⁶ DSSP assigns seven structure types (α -helix, 3(10)-helix, 5-helix, β-sheet, β -bridge, turn, and bend) according to the backbone hydrogen bonds of each residue of each conformer. If the given residue does not satisfy any of these seven types, DSSP assigns it as coil. In order to concisely summarize the differences between the force fields, we broadly classify DSSP structure types into three groups in Figure 2: (i) helix: α-helix, 3(10)-helix and 5-helix, (ii) beta: β-sheet and (iii) disordered: everything else. DSSP analysis with all eight types is presented in Figure S7–S9 for all of the proteins. Fig. 2, top row shows the average of the total structure content and it indicates that all the proteins studied are predominantly disordered with some helical content. While the helix fraction is essentially zero for CSP, it reaches 0.3 for IAPP and TDP-43_310−350. Helix fractions from the three force fields are similar to each

other, with the differences being comparable to the statistical error in the simulations. This is not unexpected, because in each of 03*, 03w and 03ws, a dihedral correction term on the ψ dihedral was optimized to reproduce helical propensity in model peptides; much larger differences would be expected when comparing to earlier force fields.^53,87 The per-residue analysis of the structure content (three bottom rows of Fig. 2) reveals that the difference in helical fractions of IAPP for different versions of the force field mainly comes from fluctuations in the middle region of the peptide, whereas the helix fractions of the two termini are mainly conserved. The TDP-43_310−350 fragment has the largest secondary structure content of all the three sequences, a significant fraction of α-helix, which is nearly identical for all three force fields (Fig. 2). CSP, on the other hand, barely has any α helix or β-sheet fraction, which is quantitatively consistent in all three force fields. In summary, therefore, three force fields with very different degrees of overall chain compaction nevertheless populate very similar secondary structures.

We note that each of these force fields have already been extensively validated against experimental NMR data for different systems. To give an illustrative example, we present NMR secondary chemical shifts and J-couplings for the TDP-43_310−350 peptide for which solution experimental data is available (Fig. 3). All calculated chemical shift deviations with respect to random coil (ΔδC_A, ΔδC_B, Δδ C, Δδ N, ΔδH_N, Δδ H_A) can be found in Figure S10,together with experimental data. We note the good consistency between the experimental and simulation CSDs (Fig. 3A) data for all force fields as evident from the average root- mean-square difference (RMSD) between experimental and simulation data for C_A− C_B (reported on legends, Fig. 3A). All of the force fields yield an RMSD (i.e. 0.94, 0.99 and 0.72 ppm), less than the prediction error of SPARTA+ (±1 ppm for carbon chemical shifts⁸⁸).CSD estimates from ff03ws are closest to the experimental data with the lowest RMSD of 0.72 ppm. Fig. 3B shows the per-residue deviation of simulated ³J_HNH coupling constant from experiment calculated for each of the three force fields, also indicating a reasonable agreement between simulation and experiment for all the three force fields. Simulated ³J_HNHA values are predicted employing the Karplus equation with a parameter set suitable for IDPs⁸⁹ (see also SI Text, Analysis). The x² difference between simulation and experiment is calculated using the experimental uncertainties,⁸⁹ the smaller the value the better the agreement between simulation and experiment. The x² from the three force fields are very similar to each other, suggesting similar performance of these force fields in reproducing the 3J_HNH coupling constant closely relevant to local torsional angles. Raw data for the three force fields and the experiment are shown in Figure S11.

Here, we emphasize that differences between the three Amber ff03 force fields (as well as the difference between the experiment and simulation) are small for both for both local structure content and NMR measurable which reflect mainly such local structure. We conclude that it is possible to match these experimental properties, essentially irrespective of how favorable the protein solvation is.

Conclusions

An accurate atomistic simulation ensemble (rigorously validated against experimental data) provides the most detailed microscopic information on IDP structure, which is diff cult to obtain solely from experimental data. We have previously shown that recent force field modifications yielded improved results in specific cases; here we have generalized to a large group of disordered proteins, finding that with the most recent ff03ws force field we are able to capture the same polymer scaling properties in water and chemical denaturant as seen in experiment. In particular, scaling exponents in water are close to 0.5, while those in high concentrations of denaturant are closer to 0.6. In contrast, older force fields yield scaling exponents in the range 0.3–0.4, characteristic of collapsed chains.

Of note is that the effects of the different force field changes appear to befforthogonal. The degree of collapse is essentially only affected by improving the solvation of the chain, here accomplished by introducing the scaling of protein-water Lennard-Jones interactions. On the other hand, however favorable the solvation model, it is still possible to reproduce local structure by altering backbone dihedral terms. Thus, dihedral angle modifications and protein solvation affect (primarily) the local structure formation and degree of collapse respectively. This stresses the need to use data reflecting both local and global structure formation in force field parameterization and validation; showing agreement with one type of data holds no guarantee for the other.

Supplementary Material