Optimized protein-water interactions and torsional refinements yield balanced atomistic protein force fields

Domain: physics
Published: 2025-11-26
Editorial outcome: major_revision

Abstract
All-atom molecular dynamics (MD) simulations based on physics-based force fields serve as an essential complement to experiments for investigating protein structure, dynamics, and interactions. Despite significant advances in force field parameterization, achieving a consistent balance of molecular interactions that stabilize folded proteins while accurately capturing the conformational dynamics of intrinsically disordered polypeptides in solution remains challenging. In this work, we introduce two refined force fields which incorporate either a selective upscaling of protein-water interactions or targeted improvements to backbone torsional sampling: (i) amber ff03w-sc, and (ii) amber ff99SBws-STQ′. Extensive validation against small-angle X-ray scattering (SAXS) and nuclear magnetic resonance (NMR) spectroscopy observables revealed that both force fields accurately reproduced the chain dimensions and secondary structure propensities of IDPs. Importantly, both force fields also maintained the stability of single-chain folded proteins and protein-protein complexes over microsecond-timescale simulations. Overall, our refinement strategies result in transferable force fields with improved accuracy for simulating diverse protein systems, ranging from folded domains to IDPs and protein-protein complexes.

Introduction
Atomistic molecular dynamics (MD) simulations based on physics-based energy functions referred to as “force fields” are now an invaluable and complementary tool to experiments in the investigation of protein structure, dynamics, and interactions associated with cellular processes1,2. With significant advances in simulation hardware and the widespread availability of conformational sampling algorithms, accessing functionally relevant biomolecular dynamics and interactions that occur on the multi-microsecond to sub-millisecond timescale are now routinely feasible3. Importantly, these advances have also enabled a rigorous evaluation of force field accuracy against experimental measurements with regard to their propensity to populate various secondary structure classes and the strength of non-bonded interactions between various chemical groups4.

After nearly two decades of force field parameterization against crystallography, spectroscopic, and ab initio calculations, modern force fields belonging to the four major families—AMBER, CHARMM, OPLS, and GROMOS—succeeded in providing a reasonable description of the structure and dynamics of folded proteins but performed poorly for short peptide ensembles in solution5,6. At present, a major emphasis and goal of modern atomistic force field development is the parameterization of transferable models which can simultaneously describe the structural stability of folded domains while capturing the transient secondary structure and global chain dimensions of intrinsically disordered polypeptides7–11. A major step towards arriving at these goals involved the extensive comparison of simulation ensembles against nuclear magnetic resonance (NMR) spectroscopy observables (chemical shifts and scalar couplings) for weakly structured peptides in solution, first carried out by Best and Hummer12. Based on these observations, global empirical corrections were applied to either ϕ or ψ backbone torsional potential for non-glycine/proline residues to adjust the helix-coil and PPII-β conformational equilibria, which led to amber ff99SB*/ff03*/ff99SBϕ13,14 and charmm22*/3615,16 force fields, which achieved an improved balance between the secondary structure classes and successfully folded both α-helical and β-sheet proteins. Notably, a similar strategy based on comparison with scalar couplings resulted in accurate sampling of side chain rotamers for isoleucine (I), leucine (L), aspartate (D), and asparagine (N) residues in ff99SB, giving rise to ff99SB*-ILDN5,15.

Despite extensive reparameterization of backbone and side chain torsion parameters, the widespread use of primitive three-site water models (e.g., TIP3P, SPC/E) consistently led to weak temperature-dependent cooperativity for protein folding17, overly collapsed structural ensembles for IDPs6,18,19, and excessive protein–protein association across the major force field families20. To remedy issues arising from weak protein–water interactions, subsequent efforts involved pairing protein force fields with accurate four-site (rigid) water models21–24. The adoption of these improved water models was a key step toward rebalancing protein–water interactions and led to improved modeling of IDP ensembles. Among these efforts, the ff99SB*/ff03* force fields were combined with the TIP4P2005 water model, along with additional upscaling (10%) of protein–water interactions and a readjustment of the global correction for ψ backbone torsion angle applied to maintain the helix-coil equilibrium for model peptides25. These modifications resulted in the ff99SBws/ff03ws force fields, which showed considerable improvements in global and structural properties computed for disordered polypeptides (from smFRET and SAXS) and reduced excessive protein–protein association25–27. These efforts were soon followed by a reparameterization of charmm36 (charmm36m) to alleviate its tendency to form left-handed α-helices, including the adoption of a modified TIP3P water with additional LJ parameters on its hydrogen atoms to enhance protein–water interactions and improve the conformational description of IDPs28. Robustelli et al. paired ff99SB-ILDN with a modified TIP4P-D water model and a Lennard-Jones (LJ) parameter modification to increase the strength of backbone hydrogen bonding, resulting in ff99SB-disp, which exhibited state-of-the-performance for both folded proteins and IDP ensembles across a wider range of test systems chosen in the study29. Along similar lines, a version of ff99SB optimized to reproduce the residue-specific torsional angle distributions observed for the PDB coil library, RSFF230, yielded highly accurate thermodynamics of protein folding when coupled with TIP4P-D water (RSFF2+)31. Other popular variants of ff99SB, such as ff14SB/19SB32,33, also alleviate excessive protein–protein association and yield more accurate IDP conformational ensembles when paired with four-site water models such as TIP4P200534 and OPC35. As an alternative to the adoption of computationally expensive four-site water models, Yoo and Aksementiev performed a reparameterization of atom pair-specific LJ parameters in charmm/amber force fields against osmotic pressure data to mitigate the over-stabilization of salt bridges and hydrophobic interactions in TIP3P water36–38.

Despite steady progress toward the development of “balanced” force fields over the last decade, numerous independent investigations reported discrepancies pertaining to certain IDP sequences, folded protein stability, and protein–protein association. For instance, ff03ws significantly overestimated (>16%) the chain dimensions of the RS peptide compared to experiment18. Further, replica-exchange simulations of WW domains indicated that ff03ws strongly destabilized the folded state37. Recent all-atom MD studies of Fused in Sarcoma (FUS), a multidomain RNA-binding protein (526 residues) with long disordered regions, also captured the structural instability of its folded RNA recognition motif (RRM) domain over the microsecond simulations with ff03ws39,40. Two independent studies that investigated the suitability of charmm36m, ff19SB-OPC, and ff99SB-disp revealed discrepancies pertaining to the solubility of Aβ peptides41 and ubiquitin self-association42, for all three force fields. These studies revealed that ff99SB-disp overestimates the strength of protein–water interactions, failing to predict β-aggregation for Aβ and weak dimerization (millimolar affinity) of ubiquitin. In contrast, charmm36m correctly predicted the aggregation behavior of Aβ wild type and two mutants, while ubiquitin self-association was found to be too strong. ff19SB-OPC exhibited intermediate behavior, giving rise to only a small population of β-aggregates for Aβ but successfully predicted weak dimerization of ubiquitin through the known interacting surfaces. To address issues related to weak protein–protein association for ff99SB-disp, Piana et al. recently proposed a reparameterization of both dihedral parameters and non-bonded interactions against osmotic pressure data43. While the reparametrized force field, referred to as DES-amber, increased the stability of protein complexes, the association free energies were still underestimated compared to experiment for some systems. Collectively, the above findings indicate that the optimal balance between protein–protein and protein-solvent interactions among modern classical force fields is yet to be fully realized.

In this work, we re-evaluate the accuracy and transferability of two Amber force fields, ff03ws and ff99SBws25, in describing the structural dynamics of both folded proteins and IDPs. We introduce two refined variants: ff03w-sc, which applies selective–water scaling to improve folded protein stability while maintaining accurate IDP ensembles, and ff99SBws-STQ′, which incorporates targeted torsional refinements of glutamine (Q) to correct overestimated helicity in polyglutamine tracts. Extensive validation against NMR and SAXS data shows that both force fields accurately reproduce IDP dimensions and secondary structure propensities, while also maintaining the stability of folded proteins and protein–protein complexes. Together, these refinements yield balanced, transferable force fields for simulating diverse protein systems.

Results
Stability of folded proteins
Recent efforts to develop “balanced” force fields—achieved by modifying protein–water vander Waals interactions or re-parameterizing water models to incorporate stronger dispersion forces—have significantly improved the prediction of global chain properties for many intrinsically disordered proteins, while maintaining essential local structural features9,23,25,26,29. However, these advancements come with a potential drawback: the strengthened protein–water interactions to better represent the IDP ensemble may inadvertently compromise the conformational stability of folded proteins23,25,29,44. To evaluate this potential issue with amber ff03ws and ff99SBws force fields, we investigated the stability of two model folded proteins: Ubiquitin (PDB ID: 1D3Z)45, a highly conserved and stable protein with well-characterized conformational dynamics across the picosecond-millisecond timescale46,47, and the Nle-double mutant of the Villin headpiece (HP35) (PDB ID: 2F4K)48, an ultrafast-folding protein48,49.

Four independent simulations (2.5 μs each) performed using ff03ws revealed significant instability for both proteins. For Ubiquitin, the most populated state deviated by approximately 0.4 nm in backbone root mean square deviation (RMSD) from the X-ray crystal structure (excluding the flexible C-terminal tail, residues 71–76), with local unfolding of the α-helix observed in one replica as indicated by per-residue root mean square fluctuation (RMSF) analysis (Supplementary Fig. 1A) and visualized over the trajectory movie (Supplementary Movie 1). Similarly, simulations of Villin HP35 exhibited pronounced structural deviations from the native state with unfolding events observed after 1 μs (Fig. 1B, Supplementary Fig. 1B, Supplementary Movie 1). In contrast, the ff99SBws force field effectively maintained the structural integrity of both Ubiquitin and Villin HP35 over the simulation timescale (Figs. 1C, D and S1). For Ubiquitin, RMSD and RMSF values remained consistently low (<0.2 nm) across all four independent simulations (Fig. 1C, Supplementary Fig. 1A, Supplementary Movie 2). Villin HP35 exhibited slightly larger structural fluctuations, with minor conformational rearrangements in the C-terminus (residues 30-35); however, the overall RMSD and RMSF remained predominantly below 0.2 nm (Figs. 1D and S1B, Supplementary Movie 2). This comparative analysis reveals that the ff99SBws force field effectively maintains native conformational states of both Ubiquitin and Villin HP35, whereas ff03ws exhibits notable destabilization with significant conformational rearrangements, indicating its potential limitations for microsecond-timescale MD simulations of folded proteins.

Modifications to the protein–water scaling scheme in Amber ff03ws
In folded proteins, structural stability is maintained through a combination of hydrogen bonding involving primarily the polypeptide backbone, and hydrophobic/electrostatic interactions mediated by side chain atoms. To investigate the relative contributions of the backbone and side chains to Ubiquitin’s stability, we initially tested two modified variants of the ff03ws force field. These variants selectively scaled protein–water LJ interactions by 10% (as described by Best et al.25) for either backbone heavy atoms (backbone scaling) or side chain atoms (side-chain scaling, excluding glycine). Similar to the original ff03ws, backbone scaling resulted in increased structural deviations (Fig. 2A, Supplementary Fig. 2C, Supplementary Movie 3), indicating a decrease in conformational stability. In contrast, side-chain scaling significantly improved the overall structural integrity of the simulation ensemble with respect to the initial structure, with RMSD fluctuations between 0.2 and 0.5 nm (Fig. 2A, Supplementary Movie 3), and RMSF values remained largely below 0.2 nm except for the disordered C-terminal region.

Charged residues are crucial for the regulation of protein stability, as they can engage in electrostatic interactions with both solvent and other residues through intra-protein salt bridges50. We therefore strategically excluded all charged residues (Lys, Arg, Asp, and Glu) from our side-chain scaling approach for several reasons. First, the parent ff03 already demonstrates robust performance in modeling salt-bridge interactions, with equilibrium association constants in good agreement with experiment, irrespective of the water model51. Second, the ff03w force field, derived from ff03* paired with the optimized TIP4P/2005 water model, yields charged side-chain hydration free energies that closely approximate experimental measurements52. Third, within our side-chain scaling framework, charged residues—which comprise approximately 30% of ubiquitin—exhibited notable deviations in solvent-accessible surface area (SASA) relative to the experimental structure (Supplementary Fig. 2). Excluding charged residues from the side-chain scaling scheme led to substantially enhanced ubiquitin stability, characterized by RMSD values consistently below 0.2 nm and reduced the solvent accessibility of charged residues (Fig. 2A, Supplementary Fig. 2C, D, Supplementary Movie 3). These improvements indicate a more stable folded conformation. This final variant, which selectively scales only uncharged side chains, is referred to as ff03w-sc and results in substantially enhanced stability for Ubiquitin.

To assess the robustness of this modified force field, we applied it to the folded state of Villin HP35, which had previously unfolded in simulations using the original ff03ws. With ff03w-sc, Villin HP35 exhibited substantially reduced structural fluctuations, with backbone RMSD and RMSF values consistently remaining below 0.2 nm (Supplementary Fig. 3A, Supplementary Movie 4), indicating improved structural stability. We further tested ff03w-sc on four additional folded proteins—GB3 (PDB ID: 1P7E)53, BPTI (PDB ID: 5PTI)54, HEWL (PDB ID: 6LYZ)55, and GTT variant of the WW domain (PDB ID: 2F21, a small β-sheet)56 —known to be destabilized by ff03ws29. While the backbone RMSD distributions for these proteins were relatively broad, they remained centered around 0.2–0.3 nm (Fig. 2B), suggesting a moderate degree of flexibility with higher RMSF values primarily localized to loop regions (Supplementary Fig. 3B–F). Overall, the proteins retained native-like conformations throughout the μs-long simulations.

To provide a more rigorous assessment of structural stability, we analyzed the fraction of native inter-residue contacts, following the approach by Best et al.57, for all six proteins using both ff03w-sc and ff99SBws. As shown in Fig. 2C and Supplementary Fig. 4, both force fields maintained high native contact fractions (>90%) across all systems. Interestingly, simulations using the ff99SBws force field consistently exhibited not only lower RMSD values (<0.2 nm) but also higher median fractions of native contacts with tighter distributions (Fig. 2B, C), indicating greater stability and closer adherence to native structures. Notably, ff99SBws also employs a 10% increase in protein–water interaction strength for all protein atoms, similar to the approach used in the original ff03ws25. Secondary structure analysis further confirmed that both α-helical and β-sheet content remained stable with minimal temporal variation throughout the simulations for all proteins (Supplementary Figs. 5 and 6). These findings demonstrate that ff03w-sc substantially improves protein structural integrity over microsecond timescales compared to the original ff03ws, while ff99SBws shows more robust performance in maintaining folded protein stability across a diverse set of test systems.

To further assess force field performance on folded proteins, we investigated their ability to capture backbone dynamics on the ps-ns timescale by calculating the amide order parameters (S²) for ubiquitin and compared them to the NMR-derived values (Supplementary Fig. 7). While the refined ff03w-sc force field successfully stabilized the overall structure of ubiquitin, it exhibited the known tendency of ff03-based force fields to overestimate backbone flexibility, particularly in loop regions where simulated S² values were consistently lower than experimental measurements (Supplementary Fig. 7A). Similar challenges have also been observed for polarizable models, where explicit inclusion of electronic polarization improves local electrostatics but does not necessarily enhance the accuracy of flexible loop dynamics, underscoring the importance of balancing bonded and non-bonded interactions58. In contrast, ff99SB-based force fields demonstrated better performance, yielding S² order parameters in excellent agreement with NMR across the entire protein and indicating more accurate descriptions of both rigid and flexible regions (Supplementary Fig. 7B). These results highlight the differential performance of the two force field families in capturing local backbone dynamics.

Chain dimensions of IDPs
To evaluate the suitability of the force fields in simulating IDPs, we determined the temperature-dependent radius of gyration (Rg) for a 34-residue fragment of Cold-shock protein (CspM34) using parallel tempering in a well-tempered ensemble (PT-WTE) to enhance conformational sampling59,60. Fig. 3A presents the Rg of CspM34 as a function of temperature for ff03w-sc and ff99SBws, compared to two reference force fields, ff03w and ff03ws. The experimental Rg estimate for CspM34 at 300 K, derived from FRET measurements using a Gaussian chain model61,62, is approximately 1.6 nm. Across the temperature range, ff03ws consistently maintained extended conformations, yielding an Rg value of 1.65 ± 0.01 nm at 300 K, which is marginally closer to the experimental measurement (~1.6 nm) than the other force fields. In contrast, ff03w produced more compact structures, with an Rg of 1.35 ± 0.1 nm at 300 K and decreasing further at elevated temperatures. Both ff03w-sc and ff99SBws yielded intermediate chain dimensions, with Rg values of 1.52 ± 0.01 nm and 1.53 ± 0.01 nm at 300 K, respectively.

We further investigated the influence of side-chain scaling on the conformational dimensions of CspM34. Variants using 10% increase in side-chain scaling, both with and without charged residues, produced comparable Rg values (1.51 ± 0.01 nm and 1.52 ± 0.01 nm at 300 K, respectively), although excluding charged residues led to slightly more collapsed conformations at higher temperatures (Supplementary Fig. 7A). These results suggest that the side-chain scaling of uncharged residues does not substantially alter the chain dimensions of CspM34. We note that 10% increase in protein side chain–water interactions used in ff03w-sc represents a near-optimal balance between maintaining folded protein structural stability and preventing over-collapse of disordered proteins. To further illustrate this trade-off, we also tested a stronger 15% scaling factor. While this higher affinity can yield more realistic chain dimensions for some disordered proteins, it also leads to the destabilization of folded proteins (Supplementary Fig. 8A–C). Therefore, although 10% scaling may not always yield the closest agreement with experimental data for every individual system, it provides reliable transferability across diverse protein systems. Based on this analysis, we adopted 10% scaling as the final choice for ff03w-sc throughout this study. Notably, this optimal scaling is specific to the base force field, as applying the same 10% side-chain adjustment to ff99SBws resulted in overly collapsed conformations of CspM34 (Supplementary Fig. 8D).

To further assess the accuracy of CspM34 ensembles generated by ff03w-sc and ff99SBws, we compared their SAXS profiles to experimental data. As an initial validation system, we selected Histatin 5, a 25-residue peptide enriched in hydrophilic amino acids (Supplementary Tables 1 and 2), which are known to promote conformational disorder compared to hydrophobic residues that are typically found in stable protein cores. Using PT-WTE simulations, we computed the SAXS profiles from the 302 K replica. Figure 3B illustrates the comparison between experimental and simulated SAXS scattering curves for Histatin 5. Both force fields demonstrated excellent agreement with experimental data, within error margins. Notably, the Rg values derived from the experimental SAXS profile63 (13.9 ± 0.1 Å) closely matched the predictions from both ff03w-sc and ff99SBws ensembles (13.8 Å).

Further validation was performed using the RS peptide and the fragment S(129–146) of the NFLt sequence (Supplementary Tables 1 and 2). Previous studies reported that ff03ws overestimated the Rg of RS peptide by ~16%18. In contrast, ff03w-sc produced SAXS profiles in excellent agreement with experimental measurements (Fig. 3C), yielding a theoretical Rg value of 12.7 ± 0.1 Å compared to the experimental value of 12.6 ± 0.1 Å. For the S(129–146) fragment64, which features alternating charged blocks that potentially facilitate salt-bridge interactions, the ff99SBws ensemble predicted a SAXS profile that aligned well with experimental data (Fig. 3D). The predicted Rg value of 12.8 ± 0.1 Å agrees closely with the experimental measurement of 12.9 ± 0.1 Å. Although ff99SB-based force fields were previously reported to overestimate salt-bridge strength51, this effect appears negligible for the S(129–146) fragment. Its moderate net negative charge (−4) and the alternating pattern of acidic (E/D) and basic (K) blocks likely promote transient rather than persistent salt bridges, while the fragment’s inherent conformational flexibility potentially further prevents their over-stabilization.

Collectively, these results demonstrate that both ff03w-sc and ff99SBws achieve robust agreement with experimental Rg values and SAXS profiles, thereby establishing their reliability in modeling the chain dimensions of small IDPs.

α-helix formation in Ace-(AAQAA)3-NH2
While accurately sampling chain dimensions in disordered peptides is crucial, it is equally important to ensure that modifications to protein–water interactions do not compromise the stability of folded motifs, such as α-helices, which may be transiently populated in unfolded states and IDPs such as the TDP-43 low-complexity domain65. To assess this balance, we investigated a 15-residue helix-forming peptide—Ace-(AAQAA)3-NH2, which exhibits approximately 30% helix content at 300 K and has a well-characterized temperature-dependent helix propensity from previous NMR studies and has become an established reference system for force field validation66. Fig. 4 illustrates the temperature-dependent helix propensity computed using ff03ws, ff03w-sc, and ff99SBws force fields. ff03w-sc, incorporating scaled side chain–water interactions, demonstrated enhanced helix stability compared to ff03ws and accurately reproduced a residual helix fraction closely aligned with NMR data at 300 K. In contrast, ff99SBws consistently overestimated the α-helical content across the temperature range and predicted elevated helicity for most residues (Fig. 4, inset). The results highlight the significant improvement achieved by ff03w-sc in stabilizing α-helix formation of the (AAQAA)3 peptide, while highlighting the need for further refinements to ff99SBws in order to accurately model the residual helicity in disordered states.

Refinement of glutamine backbone torsion parameters in ff99SBws-STQ′
Recognizing the limitations of existing force fields in accurately capturing amino acid-specific conformational preferences across diverse sequence contexts, we previously introduced targeted corrections to ff99SBws, which lowered the backbone torsion potential parameter (kψ) for serine (S), threonine (T), and glutamine (Q) from 2.0 kJ/mol to 1.0 kJ/mol, thereby decreasing the bias toward helical structures67. This refined force field—ff99SBws-STQ, substantially improved the prediction of residue-level α-helical populations in all-atom simulations of low-complexity domain (LCD) fragments derived from TAR DNA-binding protein 43 (TDP-43), FUS, and the C-terminal heptad repeat domain of RNA Polymerase II (RNAPolII)67. However, simulations of the N-terminal fragment of Httex1, comprising N17 (17 amino acids), Q16 (polyglutamine), and a minimal stretch of five prolines (P5), revealed persistent discrepancies with regard to residue-level helicity in the polyQ region: ff99SBws exhibited substantial overestimation, while ff99SBws-STQ showed significant underestimation (Fig. 5A). To address this issue, we revisited the torsion correction for glutamine and tested intermediate kψ values (between 1.0 and 2.0 kJ/mol). We found that kψ = 1.50 kJ/mol effectively corrected the bias, accurately reproducing experimental helicity in the polyQ region, with a root mean square error (RMSE) of 5.1% between simulated and experimental helix fractions (Supplementary Table 3). The revised force field, designated ff99SBws-STQ′, not only resolved the helicity discrepancies in the polyQ region but also reproduced the experimental helix fraction of the (AAQAA)3 peptide across a wide temperature range, showing excellent agreement with NMR measurements at 300 K (Fig. 5B).

To further establish the reliability of ff99SBws-STQ′, we performed a comprehensive set of validation simulations on a diverse set of both folded and disordered proteins. First, we assessed its ability to maintain the stability of folded states using the same set of six well-characterized proteins tested for ff03w-sc and ff99SBws in this study (Fig. 2B, C). These simulations yielded stable structural ensembles with RMSD values typically below 0.2 nm and a high preservation of native inter-residue contacts (>90%), comparable to those obtained with ff99SBws (Supplementary Fig. 9A, B). We then benchmarked ff99SBws-STQ′ against experimental SAXS data for the three disordered sequences: Histatin 5, RS peptide, and the S(129–146) fragment. The simulated SAXS profiles closely matched experimental scattering curves (Supplementary Fig. 9C–E), with predicted Rg values in excellent agreement with experimental measurements (Supplementary Table 4). Additionally, we evaluated 3J HN-HA scalar couplings for the RS peptide using multiple Karplus parameter sets (original Karplus68, Vogeli et al.69, DFT170, DFT270, and Larson et al.71; Supplementary Fig. 9F, Supplementary Table 5). ff99SBws-STQ′ shows the closest overall agreement with experimental values, with most residues falling within experimental uncertainty (especially with DFT1/DFT2), whereas ff99SBws slightly underestimates couplings for some central serine residues. In contrast, ff03w-sc systematically underestimates couplings across most sequence positions, consistent with prior observations for ff03-based models18. Together with global-dimension benchmarks, these results demonstrate that ff99SBws-STQ′ accurately captures global dimensions across tested systems and improves local backbone sampling for the RS peptide. However, a broader assessment of local structure across additional sequences remains an important direction for future work.

Secondary structure predictions for IDPs
To assess the applicability of ff03w-sc and ff99SBws-STQ′ for capturing transient α-helix formation in IDPs, we examined prion-like LCDs, including FUS 0−43, FUS 120−163, TDP-43 310−350, and RNA Pol II 1927−1970 72–74. We strategically selected 44-residue fragments of these proteins to optimize computational efficiency while maintaining biologically relevant regions for which detailed NMR chemical shifts and secondary structure data are readily available, allowing direct validation of our simulation results. For these fragments, we employed PT-WTE simulations to ensure converged sampling of their equilibrium properties, as demonstrated by our convergence analysis (see “Methods” and Supplementary Fig. 10). All subsequent analyses were performed using the replica at 302 K. Cα and Cβ chemical shifts were predicted using SPARTA+ for all IDP ensembles, and their conformational preferences were quantitatively assessed through secondary chemical shift difference (ΔδCα − ΔδCβ), a metric that discriminates between helical (positive values) and β-sheet (negative values) regions by comparing the observed chemical shifts against random coil reference values. Figure 6A shows that both force fields successfully captured the experimentally observed helix propensity within the core helix-forming region (residues 321–330) in TDP-43 310−350, while deviations are present in other regions of the sequence. Both force fields predicted the correct helical propensities across the chosen LCD fragments, with the exception of the N-terminal region (residue 6–13) in FUS 0−43, which was overly helical (Fig. 6), a discrepancy also reported for other force fields67. The RMSE values of simulated chemical shifts (ΔδCα − ΔδCβ) relative to experimental values for both force fields were consistently below 1 ppm (Fig. 6 and Supplementary Fig. 11), which is within the prediction error of SPARTA+75, implying the reasonable accuracy of both force fields in modeling IDP structural ensembles.

In addition, we calculated secondary structure propensities from the simulation trajectories based on the DSSP algorithm76 for each frame and averaged the helical assignments per residue over time. These propensities were compared to those predicted from experimental NMR chemical shifts using the δ2d web server77, which converts chemical shifts into residue-level secondary structure populations (Fig. 6C, D). For TDP-43, both force fields showed excellent correlation with experimental observations, accurately capturing partial helicity in the alanine-rich segment. In contrast, FUS 0−43 simulations revealed distinct behaviors: ff03w-sc exhibited elevated helical populations (~40%) in the N-terminal region (residue 6–13), while ff99SBws-STQ′ simulations predicted more transient helical conformations (<15%), closely aligning with NMR observations (Fig. 5C, D, and Supplementary Fig. 12). Analysis of additional polar-rich LCDs, including FUS 120–163 and RNA Pol II 1927−1970, demonstrated that both ff03w-sc and ff99SBws-STQ′ predicted low helical populations (<20%), consistent with experimental measurements (Supplementary Fig. 11). For reference, a comparative analysis between ff99SBws-STQ′ and ff99SB-disp (Supplementary Fig. 13) revealed that ff99SB-disp overestimated the helicity in the polyQ tract of Httex1 (residue 28–32) and the N-terminal region of FUS 0−43 (residue 6–16), while underestimating helix propensity in the helix-forming region of TDP-43 310−350 (residue 321–330).

To further evaluate the force fields for systematic, residue-specific biases, we analyzed the chemical shift deviations averaged by amino acid type across all four low-complexity fragments (Supplementary Fig. 14). This analysis shows that for most residue types, both force fields yield average deviations close to zero, indicating good overall accuracy. However, some differences emerge for polar residues. The ff03w-sc force field shows notable positive deviations for serine and threonine (similar to ff03ws67), suggesting a persistent bias toward helical conformations for these residue types. These deviations are substantially reduced in the ff99SBws-STQ′ force field, which demonstrates the success of its targeted torsional refinements for serine, threonine, and glutamine and highlights its improved accuracy for modeling polar-rich disordered sequences.

Ala 5 3J couplings
To evaluate the accuracy of backbone torsional sampling for the refined force fields, we performed 1-μs MD simulations of the alanine pentapeptide (Ala)5 and computed scalar (3J) couplings from the simulation trajectories based on the Karplus equation with different parameter sets (Table 1). The predicted values compared to the experiment were estimated using the total χ² metric (excluding 3J CC), following the procedure described by Best et al.25. A χ² value of 1.0 indicates that prediction errors match the magnitude of experimental uncertainties (σi).

Similar to ff03ws25, ff03w-sc consistently yielded χ² values below 1.0 (ranging from 0.48 to 0.80) across all tested Karplus parameter sets. These favorable χ² results (see “Methods”) suggest that side-chain-specific scaling implemented in ff03w-sc does not adversely affect backbone dihedral angle sampling, at least in small peptides. In contrast, ff99SBws-STQ′ showed slightly elevated χ² values (ranging from 1.47 to 2.09) for most Karplus parameter sets except for DFT1 (χ² = 0.99), indicating a modest overestimation of scalar couplings compared to experimental data. Collectively, these Ala5 simulations provide additional quantitative assessment for both force fields. Despite their differing refinement strategies, both force fields exhibit reasonable agreement with experimental scalar coupling data, thereby supporting their ability to achieve a reliable sampling of the conformational landscape for small peptides.

Folding miniproteins
To benchmark the transferability of the refined force fields, we evaluated their ability to reproduce the folding thermodynamics of two well-studied mini-proteins: (i) Trp-cage, a 20-residue protein containing α-helical and 310-helical elements, and (ii) chignolin, a 10-residue β-hairpin. PT-WTE folding simulations (see “Methods”) were performed across a wide temperature range using both ff03w-sc and ff99SBws-STQ′ (Supplementary Movies 5 and 6). The fraction of folded structures was calculated using a distance RMSD (dRMS) threshold of 0.2 nm from the experimental structure, as defined previously25. As shown in Fig. 7 and Supplementary Fig. 15, ff03w-sc captures the qualitative cooperative melting behavior for both mini-proteins, consistent with previous observations reported for ff03ws25, while quantitative deviations from experiment are evident across the range. In contrast, ff99SBws-STQ′ significantly underestimates the folded populations of both systems, particularly near physiological temperatures, indicating destabilization of native conformations. These results reveal that while ff99SBws-STQ′ offers improved sampling of disordered peptide ensembles, ff03w-sc maintains better thermodynamic stability for small, folded proteins. This comparison highlights the complementary strengths and inherent trade-offs of the two force fields, emphasizing the importance of selecting an appropriate force field based on the system’s structural properties.

Protein–protein complexes
Modeling protein–protein interactions within molecular complexes remains a significant challenge for most current force fields, as the multifaceted balance of water-protein, intra-protein, and water-water interactions has been primarily optimized for isolated folded or disordered proteins. Previous simulations using state-of-the-art force fields (e.g., ff99SB-disp, charmm36m, and RSFF2+) have reported instability and dissociation of many protein complexes on microsecond timescales43. To evaluate the performance of our refined force fields, we assessed the ability of ff03w-sc and ff99SBws-STQ′ to maintain stable protein–protein interfaces, focusing on four diverse complexes. Along with Barnase/Barstar (PDB ID: 1BRS)78 and SGPB/OMTKY3 (PDB ID: 3SGB)79—two systems previously identified as challenging for maintaining stability of protein–protein interactions43, we included two additional systems: CD2/CD58 (PDB ID: 1QA9)80—a β-sheet–mediated interface, and colE7/Im7 (PDB ID: 7CEI)81—a primarily α-helical complex. Each system was subjected to four independent 5-μs simulations per force field, starting from the experimentally determined structure (Supplementary Movies 7–10).

As shown in Fig. 8A, both force fields maintained the structural integrity of the complexes, with backbone RMSD values generally below 0.3 nm across all systems. ff99SBws-STQ′ consistently yielded higher populations of low-RMSD states, particularly for Barnase/Barstar and CD2/CD58, indicating enhanced structural stability. Time-evolution RMSDs of individual proteins within each complex followed similar trends (Fig. S16). To further evaluate interfacial stability, we analyzed two additional metrics: (i) the fraction of native contacts at the interface (Fig. 8B, Supplementary Fig. 17), and (ii) the interfacial surface area (Fig. 8C, Supplementary Fig. 18). Both force fields preserved high fractions of native interfacial contacts, with median stability values near 95% across the simulations. The interfacial surface areas remained tightly clustered around experimental reference values (red stars in Fig. 8C), with minimal deviations and no systematic drift over the simulation timescales.

Together, these results demonstrate that ff03w-sc and ff99SBws-STQ′ successfully maintain stable protein–protein interfaces across diverse topologies. While ff99SBws-STQ′ showed slightly improved interfacial stability, both force fields exhibited robust performance in capturing the key features of protein complex stability over extended simulation times.

Methods
Simulation of folded proteins
All-atom MD simulations of Ubiquitin, GB3, BPTI, HEWL, Villin HP35 (with double nor-leucine substitutions at residues 24 and 29), protein complexes (Barnase/Barstar and SGPB/OMTKY3) were initiated in their folded states obtained from Protein Data Bank (PDB IDs: 1D3Z, 1P7E, 5PTI, 6LYZ, 2F4K, 1BRS, and 3SGB, respectively). Systems were prepared in Gromacs 2021.590 using truncated octahedron boxes with a length of 5.5 nm, except for HEWL and protein complexes, which required a larger box of 7.0 nm and 8.0 nm, respectively. After initial energy minimization in vacuum and solvation with TIP4P/2005 water22, Na+ and Cl− ions were added to achieve 150 mM salt concentration, using improved salt parameters from Lou and Roux91. Systems were equilibrated in a canonical ensemble (NVT) using a Nose-Hoover thermostat92 with a coupling constant of 1.0 ps at 300 K, followed by an isothermal-isobaric ensemble (NPT) equilibration using a Berendsen barostat93 with anisotropic coupling constant of 5.0 ps at 1 bar. The Gromacs topologies and coordinates were then converted to Amber formats using ParmEd94, with hydrogen mass repartitioning95 to 1.5 amu to enable a 4 fs timestep. Subsequent simulations were carried out using the Amber 22 MD simulation package. The initial minimization was performed using the steepest descent and conjugate gradient algorithms, with non-hydrogen atoms restrained by a 5 kcal/mol/Å² force constant. This was followed by two 5 ns NVT equilibration phases at 300 K: the first with reduced restraints (1 kcal/mol/Å²) and the second with all restraints removed. A subsequent 10 ns NPT equilibration employed a Monte Carlo barostat (with anisotropic coupling constant of 1.0 ps and a pressure of 1.0 bar) and Langevin dynamics for temperature regulation (with a friction coefficient of 1.0 ps−1). Bonds involving hydrogen atoms were constrained using the SHAKE algorithm96. Short-range interactions were truncated at 0.9 nm, and long-range electrostatics were treated using the Particle Mesh Ewald (PME) method97,98. Finally, four independent NPT production replicas were conducted for each protein, using the same parameters as those used for equilibration.

Simulations of disordered proteins
Conformational sampling of disordered proteins was performed using PT-WTE59,60, a variant of temperature replica-exchange molecular dynamics (T-REMD)99. Simulations were conducted using Gromacs 2021.5 with the Plumed 2.8 plugin100,101. For each disordered protein sequence (Supplementary Table 1), an initial random coil conformation was generated and solvated in a truncated octahedron box containing 150 mM NaCl (using parameters from Lou and Roux) and counter ions to neutralize the system. Fourteen temperature replicas, spanning 293–500 K, were employed. Long-range electrostatics were treated with the PME method with a 0.9 nm cutoff. To enhance computational efficiency, hydrogen mass repartitioning was employed, enabling a timestep of 5 fs. Bond constraints were applied using the LINCS algorithm102. Following temperature equilibration, production runs were conducted in the NVT ensemble with biasing potentials (bias factor = 48) deposited every 4 ps. These potentials enhanced fluctuations in the system’s potential energy, resulting in a well-tempered ensemble and improved conformational exchange between replicas. The initial Gaussian potential height and width were 1.5 and 195.0 kJ/mol, respectively, with exchange probabilities between adjacent replicas ranging from 30-45%. The residue helix fraction for Httex1-Q16 was taken from Mohanty et al.103, with an aggregate trajectory of ~45 μs. Convergence of secondary structure propensities was assessed by analyzing the residual helicity at 302 K for each protein fragment. Consistent helix propensities were observed regardless of whether the entire trajectory, or trajectories excluding the initial 50 or 100 ns, were analyzed (Supplementary Fig. 11). Time-series plots of DSSP secondary structure assignments revealed transient structural behaviors, including periodic formation and melting of secondary structures throughout the trajectories (Fig. 5, Supplementary Fig. 6). These dynamic transitions facilitated convergence in estimated secondary structure propensities, underscoring the robustness of the PT-WTE approach for modeling disordered proteins.

Trajectory analyses
Structural properties were analyzed using Gromacs tools to calculate backbone RMSD, RMSF, and radius of gyration (Rg). SAXS curves were generated using the CRYSOL104, and Rg values were determined from these curves using Guinier analysis in ATLAS Primus105. Chemical shift deviations were calculated using SPARTA+75 and compared to random coil reference values obtained from the Poulsen IDP server106. Secondary structure assignments from simulations were determined using gmx do_dssp107, which implements the DSSP algorithm76. The propensity for a given structure (e.g., helix fraction) for each residue was determined by averaging the per-frame secondary structure assignments over the entire trajectory. These propensities were then compared to experimental propensities derived from NMR secondary chemical shifts using δ2d software77. SASA were calculated using the gmx sasa utility, based on the algorithm of Eisenhaber et al.108. Backbone amide S² order parameters for ubiquitin were computed using the iRED method109,110 within the Python package pyDR111,112, with the trajectories divided into 1 μs blocks for analysis as recommended by Gu et al.110. Agreement between simulated and experimental observables was quantified by calculating the RMSE using ((x_sim - x_exp)²/N)½, where N is the number of residues sampled, and x_sim and x_exp represent the simulated and experimental values, respectively. The χ² metric used to quantify agreement between simulated and experimental 3J couplings for Ala5 was defined as described previously25:

χ² = 1/N * Σ [ (J_expt,i - J_calc,i)² / σ_i² ]

where J_expt,i and J_calc,i are the experimental and computed scalar couplings, respectively, and σ_i is the estimated uncertainty associated with the predictions from the Karplus equation.

Three-dimensional structures and MD trajectories were visualized and animated using the VMD software113.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Mittal and co-workers tested the performance of two Amber force fields, namely ff03ws and ff99SBws, for both folded\ndomains and IDPs. Since ff03ws left room for improvement, they introduced AMBER ff03w-sc, which selectively scales side\nchain-water interactions for uncharged residues. The refined AMBER ff03w-sc force field yielded improved conformational\nstability of folded proteins while maintaining accurate representations of IDPs, but failed to correct the discrepancies in\nNMR-derived ps-ns timescale backbone dynamics associated with flexible loops. The performance of ff99SBws was\ngenerally better and was further improved by targeted refinement of backbone torsion parameters of serine, threonine, and\nglutamine, yielding AMBER ff99SBws-STQ’. This idea goes back to their paper in ref. 60, where they first reduced the\nbackbone torsion potential parameter k_phi from 2.0 to 1.0 kJ/mol, which was here changed to a reduction from 2.0 to 1.5\nkJ/mol.\nAll force field variants underwent thorough testing based on multiple microsecond-long MD simulations, with validations\nagainst NMR and SAXS data.\nThe manuscript is well-written, easy to follow, and features clear, interpretable figures.\nI have only two minor comments:\n1) On l. 166, the authors write that “surface-exposed charged residues” were excluded from reparameterization. I guess they\nexcluded all charged residues – there can also be charged residues more in the core of a protein, providing stability through\nsalt bridges –, and not only those on the surface. Please change the wording.\n2) The data points for the experimental results in Figs. 4, 5, and 6 are difficult to see. I recommend to plot them in yellow or\nanother bright color with sufficient contrast to the red, blue, and green colors used for the different force fields.", "reviewer_2": "In their manuscript, Phan et al. have carried out a detailed examination of the shortcomings of two all-atom force fields in the\nAMBER family (ff03ws and ff99SBws). They have also proposed improvements to each of these force fields. Their\nimprovements focus on improving sampling of realistic conformations and dynamics for folded and intrinsically disordered\nproteins, as well as protein complexes. Using both MD simulations and enhanced sampling simulation approaches,\nsystematic comparisons are made to experimental data on multiple protein systems to assess their improvements. The\nauthors identify that by scaling the sidechain-water interaction strengths—in a new force field they call ff03w-sc—they\nachieve more robust stability (specifically, RMSD) in various folded proteins relative to ff03ws. Their results suggest the\nimportance of selective modifications to uncharged residues rather than to every side chain atom. This work is an important\ncontribution to the field of biomolecular simulations, which will benefit from further work investigating the important\nproblem of force field accuracy for disordered and folded proteins.\nMajor comments:\n(1) Reproducibility\nIn the methods, insufficient detail about the simulations is provided to allow for reproducibility. For example, it is not clear in\nMethods how long the simulations are for each system. The size of each system is also not provided. In some cases, the\ninformation can be found elsewhere in the manuscript (e.g. line 142), but in other cases, it is not known. Other key\ninformation (e.g. choice of cutoff) is provided for some simulations but not others. We recommend that the others follow the\nguidelines in this checklist (Reliability and reproducibility checklist for molecular dynamics simulations. Commun Biol 6, 268\n(2023). https://doi.org/10.1038/s42003-023-04653-0) and update the methods section with complete information.\n(2) Assessment of the stability of folded proteins\nIn lines 199-200, the authors highlight “the superior performance of the original ff99sbws in maintaining the structural stability\nof folded proteins without further modification.”\nThe authors have simulated two folded proteins with this force field and evaluated stability using RMSD and RMSF. This\nevidence may not be sufficient to generalize to superior performance for all folded proteins. Furthermore, one of the authors’\nmain points is that their modified force fields maintain the stability of folded proteins. It would greatly strengthen this point to\nlook at metrics beyond RMSD/RMSF such as secondary structure, and to look at the same set of folded proteins using all\nforce fields tested.\n(3) Assessment of the accuracy of IDP ensembles\nThe authors assessed the accuracy of four small IDP ensembles (CspM34, Histatin5, RS peptide, S129-146) in terms of\nchain dimensions. However, there is additional experimental data available for these systems that could also be used to\nassess the simulation accuracy. For example, it would be useful to include the T-dependent Rg of CspM34 in Fig. 3a, as in\nRef. 25, Fig. 1. There is other experimental data for Histatin5 and RS peptide that could also be used (e.g. as in Ref. 6 & 28).\nThis would be needed to support the authors’ statement that “Further, ff99SBws also yielded structural ensembles of RS\npeptide and Histatin which are in excellent agreement with experiment.” Similarly, their statement in the Discussion and\nConclusion (“Comprehensive validation against experimental SAXS profiles and NMR data for diverse systems including\nHistatin5, RS-peptide, S(129-146) fragment, and prion-like domains - substantiates the robustness and reliability of the\nproposed refinements.”) is an overstatement given that SAXS and NMR data were not looked at for all of these systems.\n(4) Assessment of the accuracy of protein complex stability\nThe authors have assessed the stability of two protein complexes in their modified force fields based on RMSD. This\ncomparison is limited and would benefit from a deeper analysis of the interactions at the interface in both cases. Given that\nonly two complexes were studied and assessed based on RMSD alone, the statement in lines 340-342 (“These results\nhighlight the capability of ff03w-sc and ff99SBws-STQ’ to maintain the stability of protein-protein reliably interaction networks\nin complex biomolecular assemblies.”) is too strong of a claim.\n(5) Limited tests of ff99SBws-STQ’\nThe authors have carried out only limited tests of this modified force field. Lines 355-361 of the conclusion imply that\nff99SBws-STQ’ has been validated against both SAXS and NMR for all the systems listed, but this comparison has not yet\nbeen done. Further tests would be needed in order to establish the reliability of this modification. Similarly, the statement in\nthe conclusion “For\nexample, ff99SBws-STQ’ has emerged as a robust choice for simulating both folded proteins and IDPs due to its effective\nbalance of rigidity and flexibility.” is not sufficiently supported by the comparisons presented in the manuscript.\nSpecific comments:\n(1) line 64 - “…accurate four-site rigid water models(21-24) which aimed to increase the strength of the protein water\ninteractions” -> This statement implies that this was the aim of all of the studies in refs. 21-24, but that is not the case for all of\nthese studies.\n(2) lines 166-168 - “Surface-exposed charged residues are crucial to protein stability as they interact strongly with the\nsolvent and maintain essential electrostatic balances.46 We strategically excluded these residues from our side-chain\nscaling approach for several reasons.” -> this implies that it is only the surface-exposed residues that are excluded, but it’s\nactually all charged residues that are excluded.\nSimilarly, a few lines above (line 61), it is indicated that all side chain atoms have scaled interactions in the modified force\nfield.\nThe changes made to the force field become clear to the reader later in the results, but it would help to clarify the changes\nfrom the beginning.\n(3) lines 201-210 - It should be noted that the Rg of CspM34 in ff03ws is closer to the experimental value than the modified\nforce field in this case.\n(4) line 216-217 - “We note that a 10% increase in protein side chain-water interactions used in ff03w-sc represents a near-\noptimal value for balancing the structural integrity of folded proteins and preventing over-collapse of disordered proteins.” ->\nBased on Fig. S4, 15% appears to be optimal/consistent with the experiment, rather than 10%. It would be useful to comment\non this difference.\n(5) Fig. 4 — It is not clear where the experimental data shown in this figure comes from. Ref. 59 cited in the paragraph above\nFig. 4 does not show this data.\n(6) line 297-298 — “We strategically selected 44- residue fragments of these proteins to optimize computational efficiency\nwhile maintaining biological relevance and employed PT-WTE simulations to achieve converged equilibrium properties.”\nWhere is convergence of the equilibrium properties shown?\n(7) lines 302-303 — “Figure 6A shows excellent agreement between simulated and experimental data in the helix-forming\nregion (residue 321–330) in TDP-43310–350” -> this statement implies that the simulation results agree with experiment for\nall data points, which does not match what is shown in the figure.\n(8) lines 306-307 — “The RMSDs of simulated chemical shifts relative to experimental measurements for both force fields\nwere less than the prediction error of SPARTA+ (± 1 ppm for carbon chemical shifts66), highlighting their reasonable\naccuracy in modelling IDP conformational ensembles.” - It would be useful to provide the data to support this statement.\n(9) lines 316-317 — “In addition, we compared secondary structure propensities calculated from simulation trajectories using\nDSSP 67” —> DSSP as described in ref. 67 does not provide secondary structure propensities. Further details on how the\nauthors did this are needed. It would also be useful to explain why the authors chose to compare helix fraction predicted\nfrom NMR chemical shifts using a the delta2d webserver rather than directly comparing to the experimental chemical shifts.\n(10) The authors present several new results in Discussion/Conclusions that would be more appropriate in the Results\nsection (Fig. S9-S10).\n(11) Line 482 — The link to the force field files does not work.\nMinor comments:\n- To enhance the clarity of the results, the authors should state explicitly which scaling factor for the side chains is being\nused at different points in the manuscript. This becomes confusing because different possible scaling factors have been\nconsidered.\n- Line 257 — Should be Figure 4, not Figure 3.\n- Line 274-278 — Reference should be provided for this sentence\n- Line 296-298 — How were the stretches of 44 residues chosen?\n- Line 322 — Should be Figure 6c,d not 5c,d\n- Fig S3 — “ff99SBwc” in y-axis labels should be “ff99SBws”.", "reviewer_2_round_2": "The authors have addressed the majority of my earlier comments in their revised manuscript. There are a few remaining\npoints to address that could clarify and strengthen their manuscript:\n(1). Lines 237-238\n“Based on this analysis, we adopted a 10% scaling as the final choice for ff03w-sc throughout this study.”\nThe authors have partially explained why they adopted this choice of scaling. However, it should be noted that this choice\ndoes not consistently result in ensembles that are in closest agreement with the experimental data (e.g. in Figs. 3A, S8\nA&D).\n(2) Lines 325-329\n“Additionally, 3JHN-HA scalar couplings calculated for the RS-peptide showed improved agreement with experiment\nrelative to ff03w-sc and ff99SBws (Figure S9F, Table S5), implying an improvement in local backbone\nconfigurational sampling. These results demonstrate the capability of ff99SBws-STQ′ to capture both\nglobal dimensions and local structural properties across a wide range of protein systems.”\nIt is unclear what the authors mean by improved agreement here, especially for ff03w-sc. The results shown in Fig. S9F (left)\nare not in agreement with the experimental data, even qualitatively. The authors acknowledge this inconsistency in their\nresponse to the reviewers, and the details provided in the response would be beneficial to add to the manuscript.\nThe details I am referring to are:\n“Across all force fields, we evaluated couplings using multiple Karplus parameter sets, including DFT1, DFT2, and\nparameters from Vogeli et al. and Larson et al., as well as the original Karplus form. The experimental values are presented\nas a gray-shaded band (mean ± SD), facilitating direct visual comparison. As seen in the figure:\n• ff99SBws-STQ′ shows the closest agreement with experiment across the entire RS-rich sequence, with values falling\nwell within experimental uncertainty for most residues, particularly when using the DFT1 or DFT2 Karplus parameters.\n• ff99SBws also performs reasonably well, though it slightly underestimates the couplings in some central serine residues.\n• ff03w-sc only captures the qualitative trend while exhibiting systematic underestimation of couplings, with larger\ndeviations across most positions.\nWe note that this behavior of ff03w-sc is consistent with previously published results for ff03w and ff03ws, both of which\nshow similar underestimation of scalar couplings for RS-rich sequences (see Ref. 18, Figure 5A). As noted in Ref. 18, this\nsuggests that the underlying backbone φ dihedral angle preferences of the ff03-based force field family deviate from\nexperiment for this peptide, a behavior that is largely independent of the specific water-protein scaling scheme. Given that\nff03w-sc inherits its backbone parameters from ff03w, the observed behavior is in line with previous work. Overall, these\nresults support the conclusion that ff99SBws-STQ′ provides improved local conformational sampling for RS-rich\nsequences, and they further validate the force field refinement against independent NMR observables.”\n“These results demonstrate the capability of ff99SBws-STQ′ to capture both global dimensions and local structural\nproperties across a wide range of protein systems.”\nThis is an overstatement (specifically, local structural properties) given that scalar couplings were compared for only the RS\npeptide and Ala5 (and the results were not consistent with the experimental data).\n(4) Lines 393-394\n“Collectively, these Ala5 simulations provide additional quantitative validation for both force fields.”\nSimilar to the comment above, it is an overstatement to refer to this as “validation” given the inconsistency with the\nexperimental data.\n(5) Lines 446-447\n“Both force fields preserved high fractions of native interfacial contacts (>95%) across all complexes, consistently\nmaintaining key interfacial interactions.”\nThis statement is inconsistent with Fig. S17. It implies that the fraction of interfacial native contacts is always greater than\n95%.\n(6) Lines 411-413\n“As shown in Figures 7 and S15, ff03w-sc reasonably reproduces the experimental melting behavior for\nBoth mini-proteins are consistent with previous observations reported for ff03ws. Folded fractions closely\ntrack the experimental curves across the temperature range.”\nThis description is an overstatement given the clear inconsistency between the experiment (black line) and simulation (red\npoints) in Fig. 7. Lines 539-544 in the discussion would also need to be modified accordingly such that the findings are not\noverstated."}, "author_response": "Response to Referees\nWe sincerely thank both reviewers for their constructive feedback. To address their comments, we have conducted extensive\nadditional simulations and analyses, including comprehensive evaluations of both force fields across folded proteins,\ndisordered systems, and scalar coupling validations. We have also undertaken substantial revisions to improve manuscript\nclarity, enhance figure presentation, and refine our methodological descriptions. These efforts significantly strengthen the\nvalidation of our force field modifications and better align the manuscript with the reviewers’ recommendations. Below, we\nprovide detailed responses to each reviewer’s specific comments.\nReviewer #1 (Remarks to the Author):\nMittal and co-workers tested the performance of two Amber force fields, namely ff03ws and ff99SBws, for both folded\ndomains and IDPs. Since ff03ws left room for improvement, they introduced AMBER ff03w-sc, which selectively scales\nside chain-water interactions for uncharged residues. The refined AMBER ff03w-sc force field yielded improved\nconformational stability of folded proteins while maintaining accurate representations of IDPs, but failed to correct the\ndiscrepancies in NMR-derived ps-ns timescale backbone dynamics associated with flexible loops. The performance of\nff99SBws was generally better and was further improved by targeted refinement of backbone torsion parameters of serine,\nthreonine, and glutamine, yielding AMBER ff99SBws-STQ’. This idea goes back to their paper in ref. 60, where they first\nreduced the backbone torsion potential parameter k_phi from 2.0 to 1.0 kJ/mol, which was here changed to a reduction from\n2.0 to 1.5 kJ/mol.\nAll force field variants underwent thorough testing based on multiple microsecond-long MD simulations, with validations\nagainst NMR and SAXS data.\nThe manuscript is well-written, easy to follow, and features clear, interpretable figures.\nAuthors’ response: We thank the reviewer for their time and comprehensive evaluation of our manuscript. We greatly\nappreciate their insightful summary and acknowledgment of our systematic refinement efforts and extensive validation\nprocedures using multi-microsecond MD simulations, NMR, and SAXS data. We are grateful for the positive feedback and\nappreciate that the reviewer found our manuscript to be well-written, easy to follow, and supported by clear figures.\nBelow, we provide detailed responses addressing each of the specific comments and suggestions raised by the reviewer.\nI have only two minor comments:\n1) On l. 166, the authors write that “surface-exposed charged residues” were excluded from reparameterization. I guess they\nexcluded all charged residues – there can also be charged residues more in the core of a protein, providing stability through\nsalt bridges, and not only those on the surface. Please change the wording.\nAuthors’ response: We appreciate the reviewer’s insightful observation regarding the wording on line 166. The reviewer\ncorrectly noted that our initial phrasing was imprecise. Indeed, we excluded all charged residues (Lys, Arg, Asp, and Glu)\nfrom the side-chain scaling scheme, not only those exposed on the protein surface. To accurately reflect our approach, we\nhave revised the manuscript text as follows:\n“Charged residues are crucial to regulation of protein stability, as they can engage in electrostatic interactions with both\nsolvent and other residues through intra-protein salt bridges.46 We therefore strategically excluded all charged residues\n(Lys, Arg, Asp, and Glu) from our side-chain scaling approach for several reasons.”\n2) The data points for the experimental results in Figs. 4, 5, and 6 are difficult to see. I recommend to plot them in yellow\nor another bright color with sufficient contrast to the red, blue, and green colors used for the different force fields.\nAuthors’ response: We appreciate the reviewer’s suggestion regarding the visibility of experimental data points in Figures\n4, 5, and 6. Following the reviewer’s recommendation, we have updated the color of the experimental data points to yellow,\nwhich significantly enhances their contrast and visibility. These revised figures have now been incorporated into the updated\nmanuscript.\nReviewer #2 (Remarks to the Author):\nIn their manuscript, Phan et al. have carried out a detailed examination of the shortcomings of two all-atom force fields in\nthe AMBER family (ff03ws and ff99SBws). They have also proposed improvements to each of these force fields. Their\nimprovements focus on improving sampling of realistic conformations and dynamics for folded and intrinsically disordered\nproteins, as well as protein complexes. Using both MD simulations and enhanced sampling simulation approaches,\nsystematic comparisons are made to experimental data on multiple protein systems to assess their improvements. The\nauthors identify that by scaling the sidechain-water interaction strengths—in a new force field they call ff03w-sc—they\nachieve more robust stability (specifically, RMSD) in various folded proteins relative to ff03ws. Their results suggest the\nimportance of selective modifications to uncharged residues rather than to every side chain atom. This work is an important\ncontribution to the field of biomolecular simulations, which will benefit from further work investigating the important\nproblem of force field accuracy for disordered and folded proteins.\nAuthors’ response: We thank the reviewer for their thoughtful evaluation and positive feedback regarding our manuscript.\nWe are pleased that the reviewer acknowledges the systematic approach we employed to identify and address limitations\nwithin the AMBER ff03ws and ff99SBws force fields. We greatly appreciate the reviewer’s recognition of our selective\nscaling strategy for sidechain-water interactions, resulting in the ff03w-sc variant, and our comprehensive validation through\nextensive MD simulations and experimental comparisons. We are encouraged by the reviewer’s view that our work\nsignificantly contributes to advancing the accuracy of biomolecular simulations of folded proteins, intrinsically disordered\nproteins, and protein complexes.\nBelow, we provide detailed responses addressing each specific comment raised by the reviewer.\nMajor comments:\n(1) Reproducibility\nIn the methods, insufficient detail about the simulations is provided to allow for reproducibility. For example, it is not clear\nin Methods how long the simulations are for each system. The size of each system is also not provided. In some cases, the\ninformation can be found elsewhere in the manuscript (e.g. line 142), but in other cases, it is not known. Other key\ninformation (e.g. choice of cutoff) is provided for some simulations but not others. We recommend that the others follow\nthe guidelines in this checklist (Reliability and reproducibility checklist for molecular dynamics simulations. Commun Biol\n6, 268 (2023). https://doi.org/10.1038/s42003-023-04653-0) and update the methods section with complete information.\nAuthors’ response: We thank the reviewer for raising this important point regarding the clarity and reproducibility of our\nsimulation protocols. In response, we have carefully reviewed the “Reliability and reproducibility checklist for molecular\ndynamics simulations” (Commun Biol 2023) and revised our manuscript accordingly.\nWe have added a new Supplementary Table S2 that provides a comprehensive summary of the simulation conditions for all\nsystems studied. This table includes key information such as the PDB ID, box dimensions, number of atoms, number of\nreplicates, and simulation length.\nTable S2. Summary of simulation conditions for all protein systems investigated in this work.\nGroup Protein PDB ID Box size (nm) Total atoms(3) # replicas Length per replica\nUbq 1D3Z 5.5 16,403 4 2.5 μs\nHEWL 6LYZ 7.0 33,944 4 2.5 μs\nFolded\nBPTI 5PTI 5.5 16,524 4 2.5 μs\nproteins\nGB3 1P7E 5.5 16,514 4 2.5 μs\nVillin 2F4K 5.5 16,602 4 2.5 μs\nGTT 2F21 5.0 12,322 4 2.5 μs\nCspM34 N/A(1) 5.0 12,319 14 250 ns PT-WTE\n(AAQAA)\nHistatin N/A 5.0 12,376 14 250 ns PT-WTE\nRS peptide N/A 5.0 12,347 14 250 ns PT-WTE\nS(124-146) N/A 5.0 12,361 14 250 ns PT-WTE\nDisordered\nproteins\nTDP-43 N/A 6.5 27,596 14 250 ns PT-WTE\n310-350\nFUS N/A 6.5 27,567 14 250 ns PT-WTE\n0-43\nFUS N/A 6.5 27,548 14 250 ns PT-WTE\n120-163\nRNAPII N/A 6.5 27,553 14 250 ns PT-WTE\n1927-1970\nH16(2) N/A 8.0 51,577 3 250 ns PT-WTE\nJ couplings Ala N/A 4.0 6,601 1 1.0 μs\nproteins\nChignolin N/A 4.0 6,538 14 250 ns PT-WTE\nBarnase/Barstar 1BRS 8.0 50,968 4 5.0 μs\nProtein SGPB/OMTKY3 3SGB 8.0 50,958 4 5.0 μs\ncomplexes\nCD2/CD58 1QA9 8.0 50,864 4 5.0 μs\ncolE7/Im7 7CEI 8.0 50,894 4 5.0 μs\n(1)Disordered protein simulations were initiated from extended coil conformations.\n(2)The H16 simulation data using the ff99SB-disp force field were obtained from Mohanty et al. (Ref. 103), which included\nthree independent 15-μs trajectory replicates.\n(3)The reported atom counts include the dummy M sites (virtual charge sites) of the TIP4P/2005 water model.\nAdditionally, we have revised the Methods section to enhance clarity and ensured that all simulation parameters, including\nthose related to enhanced sampling and conventional MD simulations, are consistently and explicitly described.\nThese changes aim to meet the standards of transparency and reproducibility recommended by the community, and we hope\nthey will facilitate readers' ability to replicate and build upon our work.\n(2) Assessment of the stability of folded proteins\nIn lines 199-200, the authors highlight “the superior performance of the original ff99sbws in maintaining the structural\nstability of folded proteins without further modification.”\nThe authors have simulated two folded proteins with this force field and evaluated stability using RMSD and RMSF. This\nevidence may not be sufficient to generalize to superior performance for all folded proteins. Furthermore, one of the authors’\nmain points is that their modified force fields maintain the stability of folded proteins. It would greatly strengthen this point\nto look at metrics beyond RMSD/RMSF such as secondary structure, and to look at the same set of folded proteins using\nall force fields tested.\nAuthors’ response: We thank the reviewer for the thoughtful comment and the opportunity to clarify this point. While we\nappreciate the concern about generalizing conclusions from limited data, we would like to clarify that the original manuscript\ndoes not base its assessment of folded protein stability solely on two systems. In fact, we initially selected Ubiquitin and\nVillin HP35 because they are known to be challenging for the original ff03ws force field to maintain in a folded state.\nFollowing the development of ff03w-sc, we extended our benchmarking to include three additional folded proteins: HEWL\n(129 residues), BPTI (58 residues), GB3 (56 residues), as presented in Figures 2B and S3 of the original manuscript. We\nhave further expanded our dataset to a total of six proteins by now including the GTT variant of the WW domain (Figures\nS3-6, 2B and C) , a small β-sheet protein known to be unstable in ff03ws simulations (Yoon and Aksimentiev, 2016).\nTogether, these six proteins represent a diverse set in terms of size (35–129 residues), structural class (α-helical, β-sheet,\nand mixed α/β), and sequence complexity. Further, this dataset includes both proteins that were previously shown to\nhighlight weaknesses of existing force fields (ff03ws and ff9SBws) and those commonly used as benchmarks for structural\nstability.\nFigure 2: Stability of folded proteins using ff03w-sc and ff99SBws force fields. (A) Time-dependent backbone RMSD of ubiquitin\nfor three test scaling schemes derived from ff03ws: backbone scaling, sidechain scaling for all non-glycine residues, and final ff03w-sc\nscheme with sidechain scaling for uncharged residues only. The red dashed line at RMSD = 0.2 nm corresponds to a stable conformation.\n(B) Distribution of backbone RMSD values for six folded proteins—Ubiquitin (Ubq), HEWL, BPTI, GB3, Villin HP35, and GTT variant\nof WW domain—simulated using ff03w-sc (blue) and ff99SBws (yellow). Distributions reflect combined data from four independent\nsimulations for each system. (C) Fraction of native contacts retained over the course of simulations, calculated using the method of Best,\nHummer, and Eaton (2013).53 Boxplots represent the distribution across four independent simulations per system.\nIn addition to RMSD and RMSF, we have now incorporated additional metrics to provide a more comprehensive assessment\nof structural stability:\n• Native contact analysis: We calculated the fraction of native contacts (as defined by Best, Hummer, and Eaton, 2013) for\nall six proteins. As shown in the revised Figure 2C, both ff03w-sc and ff99SBws maintain high native contact fractions\n(>90%) across all systems. ff99SBws shows notably better performance with a higher median fraction of native contacts\nand tighter distributions, indicating less deviation from native states.\n• Secondary structure stability: We have also analyzed the time evolution of secondary structure using DSSP assignments\nand computed secondary structure propensities over time. The results, presented in new Supporting Information Figures S5-\n6, show that α-helical and β-sheet content remains highly stable for all proteins throughout multi-microsecond simulations,\nconfirming that the structural cores of the folded states are well preserved.\nWe have revised the manuscript to incorporate this more comprehensive validation with revised text additions and changes\nas provided below:\n“...We further tested ff03w-sc on four additional folded proteins—GB3 (PDP ID: 1P7E),51 BPTI (PDB ID: 5PTI),52 and\nHEWL (PDB ID: 6LYZ),53 and GTT variant of the WW domain (PDB ID: 2F21, a small β-sheet),54—known to destabilize\nin ff03ws simulations.29….To provide a more rigorous assessment of structural stability, we analyzed the fraction of native\ninter-residue contacts, following the approach by Best et al.,53 for all six proteins using both ff03w-sc and ff99SBws . As\nshown in Figure 2C, both force fields maintained high native contact fractions (>90%) across all systems. Interestingly,\nsimulations using the ff99SBws force field consistently exhibited not only lower RMSD values (<0.2 nm) and narrower\ndistributions but also higher median fractions of native contacts with tighter distributions (Figures 2B, C), indicating\ngreater stability and closer adherence to native structures. Notably, ff99SBws also employs a 10% increase in protein-\nwater interaction strength for all protein atoms, similar to the approach used in the original ff03ws.25 Secondary structure\nanalysis further confirmed that both α-helical and β-sheet content remained stable with minimal temporal variation\nthroughout the simulations for all proteins (Figures S5-6). These findings demonstrate that ff03w-sc substantially improves\nprotein structural integrity over microsecond timescales compared to the original ff03ws, while ff99SBws shows more\nrobust performance in maintaining folded protein stability across a diverse set of test systems.”\n(3) Assessment of the accuracy of IDP ensembles\nThe authors assessed the accuracy of four small IDP ensembles (CspM34, Histatin5, RS peptide, S129-146) in terms of\nchain dimensions. However, there is additional experimental data available for these systems that could also be used to\nassess the simulation accuracy. For example, it would be useful to include the T-dependent Rg of CspM34 in Fig. 3a, as in\nRef. 25, Fig. 1. There is other experimental data for Histatin5 and RS peptide that could also be used (e.g. as in Ref. 6 &\n28). This would be needed to support the authors’ statement that “Further, ff99SBws also yielded structural ensembles of\nRS peptide and Histatin which are in excellent agreement with experiment.” Similarly, their statement in the Discussion\nand Conclusion (“Comprehensive validation against experimental SAXS profiles and NMR data for diverse systems\nincluding Histatin5, RS-peptide, S(129-146) fragment, and prion-like domains - substantiates the robustness and reliability\nof the proposed refinements.”) is an overstatement given that SAXS and NMR data were not looked at for all of these\nsystems.\nAuthors’ response: We thank the reviewer for this thoughtful and detailed comment. We appreciate the opportunity to clarify\nand improve the manuscript.\n1. Temperature-dependent R of CspM34:\ng\nWe have updated Figures 3A, S8A,D to include the experimental temperature-dependent R values of CspM34 from Wuttke\ng\net al. (Ref. 59). This analysis provides a quantitative comparison to our simulation results across the entire temperature\nrange. We have also revised the figure caption to reflect this addition. As evident from the plot provided below, adoption of\nthe sidechain scaling scheme (ff03ws-sc) does not have a significant effect on the temperature dependent R profile\ng\ncompared to the parent force field (ff03ws).\nFigure 3. Chain dimensions of IDPs. (A) Radius of gyration of CspM34 as a function of temperature. Results are shown for\nff03ws (green), ff03w-sc (blue), ff99SBws (yellow), and ff03w (black), respectively. Yellow plus signs indicate values estimated\nfrom FRET measurements.59 The experimental value at 300 K is indicated by the red dashed line. (B) …\n2. Validation of RS-peptide using ³JHN–HA coupling constants:\nIn response to the reviewer’s suggestion, we performed a detailed assessment of our simulation ensemble against ³J\nHN–HA\ncoupling constants, which are highly sensitive to the backbone φ dihedral angle distribution. We calculated the ³J\nHN-\ncouplings for the RS-peptide and compared them to experimental values from Ref. 18. The analysis is shown in the new\nHA\nSupplementary Figure S9F for three different force fields: ff03w-sc, ff99SBws, and ff99SBws-STQ′.\nFigure S9. …(F) ³J scalar couplings for the RS-peptide, computed using multiple Karplus parameter sets for simulations with\nHN–HA\nff03w-sc, ff99SBws, and ff99SBws-STQ′. Experimental values (gray band) are from Ref. 18.\nAcross all force fields, we evaluated couplings using multiple Karplus parameter sets, including DFT1, DFT2, and\nparameters from Vogeli et al. and Larson et al., as well as the original Karplus form. The experimental values are presented\nas a gray-shaded band (mean ± SD), facilitating direct visual comparison. As seen in the figure:\n• ff99SBws-STQ′ shows the closest agreement with experiment across the entire RS-rich sequence, with values falling\nwell within experimental uncertainty for most residues, particularly when using the DFT1 or DFT2 Karplus parameters.\n• ff99SBws also performs reasonably well, though it slightly underestimates the couplings in some central serine residues.\n• ff03w-sc only captures the qualitative trend while exhibiting systematic underestimation of couplings, with larger\ndeviations across most positions.\nWe note that this behavior of ff03w-sc is consistent with previously published results for ff03w and ff03ws, both of which\nshow similar underestimation of scalar couplings for RS-rich sequences (see Ref. 18, Figure 5A). As noted in Ref. 18, this\nsuggests that the underlying backbone φ dihedral angle preferences of the ff03-based force field family deviate from\nexperiment for this peptide, a behavior that is largely independent of the specific water-protein scaling scheme. Given that\nff03w-sc inherits its backbone parameters from ff03w, the observed behavior is in line with previous work. Overall, these\nresults support the conclusion that ff99SBws-STQ′ provides improved local conformational sampling for RS-rich\nsequences, and they further validate the force field refinement against independent NMR observables.\nWe have added a summary of this analysis to the Results section as shown below and provided the new figure in the\nSupporting Information.\n“Additionally, 3J scalar couplings calculated for the RS-peptide showed improved agreement with experiment relative\nHN-HA\nto ff03w-sc and ff99SBws (Figure S9F, Table S5), supporting enhanced local backbone sampling. These results demonstrate\nthe capability of ff99SBws-STQ′ to capture both global dimensions and local structural properties across a wide range of\nprotein systems.”\n3. Clarification on the summary statement - “Comprehensive validation against experimental SAXS profiles and NMR\ndata for diverse systems including Histatin5, RS-peptide, S(129-146) fragment, and prion-like domains - substantiates the\nrobustness and reliability of the proposed refinements” in the Discussion and Conclusion.\nWe appreciate the reviewer’s careful reading and agree that our original wording was an oversight in our writing. Our\nintention was to specify which experimental data corresponded to which systems. We have revised the sentence in the\nDiscussion and Conclusion to accurately reflect the comparisons performed:\n“Validation against experimental SAXS profiles for Histatin5, RS-peptide, and the S(129–146) fragment, along with NMR-\nderived secondary structure for prion-like domains and scalar coupling analysis for RS-peptide and Ala , supports the\n5\nrobustness and transferability of the proposed force field refinements.”\nAdditionally, please see our response to point 5 of the major comment regarding the scalar J-coupling analysis for Ala .\n5\n(4) Assessment of the accuracy of protein complex stability\nThe authors have assessed the stability of two protein complexes in their modified force fields based on RMSD. This\ncomparison is limited and would benefit from a deeper analysis of the interactions at the interface in both cases. Given that\nonly two complexes were studied and assessed based on RMSD alone, the statement in lines 340-342 (“These results\nhighlight the capability of ff03w-sc and ff99SBws-STQ’ to maintain the stability of protein-protein reliably interaction\nnetworks in complex biomolecular assemblies.”) is too strong of a claim.\nAuthors’ response: We thank the reviewer for this valuable suggestion to expand and deepen our assessment of protein\ncomplex stability. In response, we have performed additional simulations and interfacial analyses to provide a more\ncomprehensive evaluation.\nSpecifically, we added two additional protein complexes to complement our original study of Barnase/Barstar and\nSGPB/OMTKY3: CD2/CD58 (PDB ID 1QA9), a β–sheet–mediated interface, and colE7/Im7 (PDB ID 7CEI), a primarily\nα-helical interface. Each system was simulated using four independent 5-μs trajectories for both ff03w-sc and ff99SBws-\nSTQ’. The new results are presented in the revised Figure 7 and new supplementary figures (Figures S10-12).\nFigure 8. Protein complex stability. Analyses of four protein complexes simulated using ff03w-sc and ff99SBws-STQ' force fields.\nThe insets present the experimental structures of the complexes. (A) Distributions of the backbone RMSD of the full complexes. (B)\nAverage fraction of native contacts at the interface; error bars represent the standard deviation across four independent simulations. (C)\nViolin plots showing the distribution of the interfacial surface area; the red stars indicate the values calculated from the experimental\ncrystal structures.\nIn terms of backbone RMSD distributions of the overall complexes (Figure 8A), ff99SBws-STQ ensembles consistently\nexhibit a higher population of low-RMSD states, particularly for the β-sheet-rich CD2/CD58 and Barnase/Barstar\ncomplexes, suggesting slightly enhanced structural stability at the interface. ff03w-sc also preserved overall complex\nintegrity across all systems, though with slightly broader RMSD distributions, generally remaining below 0.3–0.35 nm.\nTime-dependent RMSDs of individual proteins within each complex showed similar trends (Figure S16).\nIn addition to RMSD, we analyzed two critical interfacial stability metrics:\n● Fraction of native contacts (Figures 7B, S17), which reflects the preservation of specific residue-level interactions at\nthe protein–protein interface.\n● Interfacial surface area (Figures 7C, S18), which reports the extent of buried surface area upon complex formation.\nAs shown in Figure 8B, both force fields successfully maintained a high fraction of native contacts (>95%) across all four\ncomplexes, indicating consistent preservation of key interfacial interactions. Furthermore, the interfacial surface areas\nremained tightly clustered around the experimental reference values (Figure 8C), with only minor deviations and no\nsystematic drift. These findings indicate that both ff03w-sc and ff99SBws-STQ can maintain stable interface geometries\nover long simulation timescales.\nWe have updated the “Protein–Protein Complexes” section of the manuscript to incorporate these new results and provide\na more balanced and comprehensive evaluation of complex stability as provided below:\n“Protein–Protein Complexes. Modeling protein–protein interactions within molecular complexes remains a significant\nchallenge for many current force fields, as the balance among water–protein, intra-protein, and water–water interactions is\noften optimized for isolated or disordered proteins. Prior studies using state-of-the-art force fields (e.g., ff99SB-disp,\nCHARMM36m, and RSFF2+) have reported partial unfolding or dissociation of several protein complexes on the\nmicrosecond timescale.41 To evaluate the performance of our refined force fields, we assessed the ability of ff03w-sc and\nff99SBws-STQ to maintain stable protein–protein interfaces, focusing on four diverse complexes. Along with\nBarnase/Barstar (PDB ID: 1BRS)79 and SGPB/OMTKY3 (PDB ID: 3SGB)80 - two systems previously identified as\nchallenging for maintaining stability of protein-protein interactions,41 we included two additional systems: CD2/CD58\n(PDB ID: 1QA9)81 - a β-sheet–mediated interface, and colE7/Im7 (PDB ID: 7CEI)82 - a primarily α-helical complex. Each\nsystem was subjected to four independent 5-μs simulations per force field, starting from the experimentally determined\nstructure (Movies S7-10).\nAs shown in Figure 8A, both force fields maintained the structural integrity of the complexes, with backbone RMSD\nvalues generally below 0.3 nm across all systems. ff99SBws-STQ consistently yielded higher populations of low-RMSD\nstates, particularly for Barnase/Barstar and CD2/CD58, indicating enhanced structural stability. Time-evolution RMSDs\nof individual proteins within each complex followed similar trends (Fig. S16). To further evaluate interfacial stability, we\nanalyzed two additional metrics: (i) the fraction of native contacts at the interface (Figures 8B, S17), and (ii) the interfacial\nsurface area (Figures 8C, S18). Both force fields preserved high fractions of native interfacial contacts (>95%) across all\ncomplexes, consistently maintaining key interfacial interactions. The interfacial surface areas remained tightly clustered\naround experimental reference values (red stars in Figure 8C), with minimal deviations and no systematic drift over the\nsimulation timescales.\nTogether, these results demonstrate that ff03w-sc and ff99SBws-STQ successfully maintain stable protein–protein\ninterfaces across diverse topologies. While ff99SBws-STQ showed slightly improved interfacial stability, both force fields\nexhibited robust performance in capturing the key features of protein complex stability over extended simulation times.”\n(5) Limited tests of ff99SBws-STQ’\nThe authors have carried out only limited tests of this modified force field. Lines 355-361 of the conclusion imply that\nff99SBws-STQ’ has been validated against both SAXS and NMR for all the systems listed, but this comparison has not yet\nbeen done. Further tests would be needed in order to establish the reliability of this modification. Similarly, the statement\nin the conclusion “For example, ff99SBws-STQ’ has emerged as a robust choice for simulating both folded proteins and\nIDPs due to its effective balance of rigidity and flexibility.” is not sufficiently supported by the comparisons presented in\nthe manuscript.\nAuthors’ response: We thank the reviewer for this insightful comment. To address this concern, we have conducted a\ncomprehensive set of new simulations and analyses for both folded and disordered proteins using the ff99SBws-STQ′ and\nff03ws-sc force fields.\n1. Validation on folded proteins: We simulated a diverse set of six well-folded proteins—Ubiquitin, HEWL, BPTI, GB3,\nVillin HP35, and the GTT variant of the WW domain—using the ff99SBws-STQ force field. This dataset spans different\nsizes (35–129 residues) and structural motifs (α-helical, β-sheet, and mixed topologies). As shown in Figures S9A–B, the\nRMSD distributions remain tightly centered below 0.2–0.25 nm, and the fraction of native contacts is consistently high\n(>90%) across all proteins. These metrics indicate that ff99SBws-STQ′ maintains the integrity of the folded state\ncomparably to ff99SBws.\n2. Validation on disordered proteins: We also simulated three disordered sequences—Histatin5, RS-peptide, and the\nS(129–146) fragment—and compared the predicted SAXS profiles (global dimensions) with experimental data. As shown\nin Figures S9C–E, the simulated scattering curves (red dashed lines) closely match with experiment (gray), and the predicted\nradii of gyration (13.8 ± 0.1 Å for Histatin5, 12.7 ± 0.1 Å for RS-peptide, and 12.5 ± 0.1 Å for S129–146) are in good\nagreement with reported experimental values (13.9 ± 0.1, 12.6 ± 0.1, and 12.9 ± 0.1 Å, respectively).\nIn addition, we computed 3J scalar couplings for the RS-peptide based on multiple Karplus parameter sets. As shown\nHN–HA\nin Figure S9F, ff99SBws-STQ′ predicted coupling constants that closely match experimental values across the sequence,\nwith improved accuracy relative to ff03w-sc and ff99SBws. This demonstrates improved sampling of backbone φ dihedral\nangles. For additional details, please refer to our response to point 3 of the major comments.\nFigure S9. Validation of folded and disordered protein simulations using ff99SBws-STQ′. (A) Backbone RMSD distributions for\nsix folded proteins—Ubiquitin (Ubq), HEWL, BPTI, GB3, Villin HP35, and the GTT variant of the WW domain—simulated using\nff99SBws (yellow) and ff99SBws-STQ′ (red). RMSD values are calculated relative to the experimental structures, with data combined\nfrom four independent simulations for each system. (B) Fraction of native inter-residue contacts retained throughout the simulations for\nthe same six proteins. Boxplots summarize the distribution across simulation replicates. (C–E) Comparison of simulated SAXS profiles\n(red dashed lines) and radii of gyration to experimental data (gray points with error bars) for three disordered proteins—Histatin5, RS-\npeptide, and S(129–146)—simulated with ff99SBws-STQ′. Predicted Rg values (13.8 ± 0.1 nm for Histatin5, 12.7 ± 0.1 nm for RS-\npeptide, and 12.5 ± 0.1 nm for S(129–146)) are in close agreement with experimental values (13.9 ± 0.1, 12.6 ± 0.1, and 12.9 ± 0.1 nm,\nrespectively). (F) ³JHN–HA scalar couplings for the RS-peptide, computed using multiple Karplus parameter sets for simulations with\nff03w-sc, ff99SBws, and ff99SBws-STQ′. Experimental values (gray band) are from Ref. 18.\nWe have added a new paragraph in the Results section to incorporate these results.\n“To further establish the reliability of ff99SBws-STQ′ , we performed a comprehensive set of validation simulations on a\ndiverse set of both folded and disordered proteins. First, we assessed its ability to maintain the stability of folded states\nusing the same set of six well-characterized proteins tested for ff03w-sc and ff99SBws in this study (Figure 2B, C). These\nsimulations yielded stable structural ensembles with RMSD values typically below 0.2 nm and a high preservation of native\ninter-residue contacts (>90%), comparable to those obtained with ff99SBws (Figures S9A-B). We then benchmarked\nff99SBws-STQ′ against experimental SAXS data for the three disordered sequences: Histatin5, RS-peptide, and the S(129–\n146) fragment. The simulated SAXS profiles closely matched experimental scattering curves (Figure S9C-E), with predicted\nR values in excellent agreement with experimental measurements (Table S4). Additionally, 3J scalar couplings\ng HN-HA\ncalculated for the RS-peptide showed improved agreement with experiment relative to ff03w-sc and ff99SBws (Figure S9F,\nTable S5), implying an improvement in local backbone configurational sampling. These results demonstrate the capability\nof ff99SBws-STQ′ to capture both global dimensions and local structural properties across a wide range of protein systems.”\n3. Validation on the Ala5 system. To provide additional quantitative validation, we also performed 1-μs simulations of the\nAlanine pentapeptide (Ala5), a well-characterized model system and evaluated the agreement between simulated and\nexperimental scalar couplings using the total χ² metric (excluding ³J ), following the procedure described by Best et al.\nCC\n(2014). For a perfect force field, χ² value of 1.0 is expected which indicates that the prediction error for the Karplus equation\nis equal to the experimental uncertainty (σ).\ni\nSimilar to the parent force field (Best et al., 2014), ff03w-sc demonstrates χ² values less than 1.0 across all parameter sets.\nAs suggested previously, χ² < 1 likely arises due to an overestimation of the experimental σ which are obtained from diverse\ni\nsequences in addition to olgio-alanine. Since ff03w-sc retains the original ff03w backbone potential, these favorable χ²\nresults suggest that its side-chain–specific scaling does not adversely affect backbone sampling, at least in the context of\nsmall peptides. With the exception of DFT1, ff99SBws-STQ′ yields χ² values > 1 (ranging from 1.47 to 2.09) across all\nother Karplus parameter sets, which implies a modest overestimation of the coupling constants compared to experiment.\nTogether, these Ala simulations provide additional validation that both force fields—despite differing strategies for\n5\nrefinement—achieve reasonable agreement with experimental scalar coupling data and further highlight their ability to\nreliably capture the observed conformational behavior for small peptides. We have added a paragraph on “Ala 3J couplings”\n5\nin the Results section.\n“Ala 3J couplings. To evaluate the accuracy of backbone torsional sampling for the refined force fields, we performed 1-\n5\nμs molecular dynamics simulations of the alanine pentapeptide (Ala ) and computed \\′. scalar (3J) couplings from the\n5\nsimulation trajectories based on the Karplus equation with different parameter sets. The predicted values were compared\nto experiment were estimated using the total χ² metric (excluding 3J ), following the procedure described by Best et al.\nCC\n(2014).25 A χ² value of 1.0 indicates that prediction errors match the magnitude of experimental uncertainties (σ).\ni\nTable 1. Total χ² values† of 3J couplings (excluding 3JCC) of Ala .\nOriginal73 0.73 1.68\nVogeli et al.74 0.54 2.09\nDFT175 0.80 0.99\nDFT275 0.48 1.47\nLarson et al.76 0.77 1.72\n†A full list of 3J couplings for each residue is given in Table S6.\nSimilar to ff03ws,25 ff03w-sc consistently yielded χ² values below 1.0 (ranging from 0.48 to 0.80) across all tested\nKarplus parameter sets. These favorable χ² results (see Methods) suggest that side-chain–specific scaling implemented in\nff03w-sc does not adversely affect backbone dihedral angle sampling, at least in small peptides. In contrast, ff99SBws-STQ′\nshowed slightly elevated χ² values (ranging from 1.47 to 2.09) for most Karplus parameter sets except for DFT1 (χ² = 0.99),\nindicating a modest overestimation of scalar couplings compared to experimental data. Collectively, these Ala5 simulations\nprovide additional quantitative validation for both force fields. Despite their differing refinement strategies, both force\nfields exhibit reasonable agreement with experimental scalar coupling data, thereby confirming their ability to achieve a\nreliable sampling of the conformational landscape for small peptides.”\n4. Validation on protein folding thermodynamics. Finally, to evaluate the performance of the refined force fields on the\ncomplex process of protein folding , we investigated the folding equilibria of two canonical mini-proteins: (i) the 20-residue\nTrp-cage, which adopts an α-helical and 3 -helical structure, and (ii) the 10-residue chignolin, which forms a β-hairpin\n10\nstabilized by turns and backbone hydrogen bonding. Replica-exchange folding simulations were conducted across a wide\ntemperature range using ff03w-sc and ff99SBws-STQ′, and the folded fraction was computed based on a distance RMSD\n(dRMS) threshold of 0.2 nm from the experimental structure, following the definition in Best et al. (2014).\nAs shown in Figures 7 and S15, ff03w-sc (blue squares) recapitulates the melting behavior of both Trp-cage and chignolin,\nwith folded fractions closely tracking experimental melting curves (black lines) across the temperature range. In contrast,\nff99SBws-STQ′ (red circles) significantly underestimates the folded population at all temperatures for both systems,\nsuggesting that it may overly destabilize compact helical and β-structured conformations in small proteins. These results\nhighlight a key distinction between the force fields: while ff99SBws-STQ′ improves local sampling and disordered protein\nbehavior through torsional refinements, ff03w-sc offers better thermodynamic stability for small, folded proteins.\nThese folding benchmarks offer additional insight into the strengths and trade-offs of each refined force field, highlighting\nthe importance of selecting an appropriate model based on the dominant structural features of the system under study.\nWe have incorporated these findings into a new paragraph, “Folding mini-proteins”, in the Results section.\n“Folding mini-proteins. To benchmark the transferability of the refined force fields, we evaluated their ability to reproduce\nthe folding thermodynamics of two well-studied mini-proteins: (i) Trp-cage, a 20-residue protein containing α-helical and\n3 -helical elements, and (ii) chignolin, a 10-residue β-hairpin. PT-WTE folding simulations (see Methods) were performed\n10\nacross a wide temperature range using both ff03w-sc and ff99SBws-STQ′ (Movies 5-6). The fraction of folded structures\nwas calculated using a distance RMSD (dRMS) threshold of 0.2 nm from the experimental structure, as defined previously.25\nAs shown in Figures 7 and S15, ff03w-sc reasonably reproduces the experimental melting behavior for both mini-\nproteins, consistent with previous observations reported for ff03ws.25 Folded fractions closely track the experimental curves\nacross the temperature range. In contrast, ff99SBws-STQ′ significantly underestimates the folded populations of both\nsystems, particularly near physiological temperatures, indicating destabilization of native conformations. These results\nreveal that while ff99SBws-STQ′ offers improved sampling of disordered peptide ensembles, ff03w-sc maintains better\nthermodynamic stability for small, folded proteins, likely due to the removal of backbone-water interaction scaling in\ncombination with the original ff03w backbone torsional modification applied to reproduce the helix-coil equilibria for\nshort peptides in TIP4P/2005 water. This comparison highlights the complementary strengths and inherent trade-offs of the\ntwo force fields, emphasizing the importance of selecting an appropriate force field based on the system’s structural\nproperties.\nFigure 7. Calculated melting curves of mini-proteins simulated with ff03w-sc and ff99SBws-STQ′. (A) Temperature-dependent folded\nfraction of Trp-cage and (B) chignolin, computed from replica exchange simulations using ff03w-sc (blue squares) and ff99SBws-STQ′\n(red circles). A structure is classified as folded if its distance RMSD (dRMS) from the experimental structure is less than 0.2 nm, as\ndefined in Best et al. (2014). Error bars represent standard deviations obtained from block averaging over five equal-length blocks. The\nexperimental melting curves (black lines) are included for comparison. Insets show the representative folded structures of Trp-cage\n(PDB ID: 1L2Y)77 and chignolin (PDB ID: 1UAO)78.\nWe have also revised the “Discussion and Conclusion” section to incorporate these new results and provide a more balanced\nassessment of force field performance. Importantly, both ff03w-sc and ff99SBws-STQ′ have now been systematically\nevaluated across the complete range of protein systems presented in this study—including folded proteins (six well-\ncharacterized systems), disordered proteins (three sequences with SAXS validation and four 44-residue LCDs with NMR\nvalidation), mini-protein folding (Trp-cage and chignolin), small peptide conformational sampling (Ala scalar couplings),\n5\nand protein-protein complexes (four systems with α-helical-mediated, β-sheet-mediated, and mixed interfaces)—thereby\nestablishing a comprehensive benchmark that provides robust validation for both force field refinements.\nSpecific comments:\n(1) line 64 - “…accurate four-site rigid water models(21-24) which aimed to increase the strength of the protein water\ninteractions” -> This statement implies that this was the aim of all of the studies in refs. 21-24, but that is not the case for\nall of these studies.\nAuthors’ response: We thank the reviewer for this important correction. The reviewer is correct that our original wording\nwas an overgeneralization, as not all studies cited in Refs. 21–24 were specifically aimed at increasing the strength of\nprotein-water interactions. Our intention was to highlight that the adoption of these more accurate water models was a key\nstrategy used by the field to improve the description of protein solvation. We have revised the sentence for clarity.\n“To remedy issues arising from weak protein-water interactions, subsequent efforts involved pairing protein force fields\nwith accurate four-site (rigid) water models.21–24 The adoption of these improved water models was a key step toward\nrebalancing protein-water interactions and led to improved modeling of IDP ensembles.”\n(2) lines 166-168 - “Surface-exposed charged residues are crucial to protein stability as they interact strongly with the\nsolvent and maintain essential electrostatic balances.46 We strategically excluded these residues from our side-chain scaling\napproach for several reasons.” -> this implies that it is only the surface-exposed residues that are excluded, but it’s actually\nall charged residues that are excluded.\nSimilarly, a few lines above (line 61), it is indicated that all side chain atoms have scaled interactions in the modified force\nfield.\nThe changes made to the force field become clear to the reader later in the results, but it would help to clarify the changes\nfrom the beginning.\nAuthors’ response: We appreciate the reviewer’s comment on improving the clarity of our manuscript. As noted in our\nresponse to point 1 from Reviewer 1, we have revised the relevant sentence to clarify that all charged residues (Lys, Arg,\nAsp, and Glu) were excluded from the side-chain scaling scheme, not only those exposed on the protein surface.\nRegarding the second point, the reviewer is correct that our initial description suggested all side chains were scaled but later\nclarifies that only uncharged residues were included in the final implementation. This progression was intentional as we\naimed to present the logical sequence of modifications that led to the development of the ff03w-sc force field. To improve\nclarity while preserving this stepwise narrative, we have revised the relevant text in the Results sections to better guide the\nreader through the evolution of the force field and explicitly distinguish exploratory tests from the final model.\nParagraph starts at line 156:\n“…To investigate the relative contributions of the backbone and side chain interactions towards Ubiquitin’s stability, we\ninitially tested two modified variants of the ff03ws force field…”\nParagraph starts at line 166:\n“...Excluding charged residues from the side-chain scaling scheme led to substantially enhanced ubiquitin stability,\ncharacterized by RMSD values consistently below 0.2 nm and reduced the solvent accessibility of charged residues (Figures\n2A, S2C, D, Movie S3). This final variant, which selectively scales only uncharged side chains, is referred to as ff03w-sc\nand results in substantially enhanced stability for Ubiquitin.”\n(3) lines 201-210 - It should be noted that the Rg of CspM34 in ff03ws is closer to the experimental value than the\nmodified force field in this case.\nAuthors’ response: We thank the reviewer for this careful observation regarding the radius of gyration (R ) of CspM34. The\ng\nreviewer is correct that, for this metric, the original ff03ws force field yields an R value (1.65 ± 0.01 nm) that is marginally\ng\ncloser to the experimental measurement (~1.6 nm) than the value obtained with the modified ff03w-sc force field (1.52 ±\n0.01 nm). However, the primary motivation for developing ff03w-sc was to overcome a critical limitation of ff03ws: its\ninability to maintain the structural stability of folded proteins. While ff03ws performs reasonably well for disordered\nproteins, it exhibits significant destabilization in simulations of folded proteins. In contrast, ff03w-sc offers a more balanced\nperformance across both folded and disordered systems.\nTo improve clarity, we have revised the relevant sentence to note this comparison in the chain dimensions of CspM34.\n“Across the temperature range, ff03ws consistently maintained extended conformations, with an Rg value of 1.65 ± 0.01\nnm at 300 K, which is marginally closer to the experimental measurement (~1.6 nm) than the other force fields.”\n(4) line 216-217 - “We note that a 10% increase in protein side chain-water interactions used in ff03w-sc represents a near-\noptimal value for balancing the structural integrity of folded proteins and preventing over-collapse of disordered proteins.”\n-> Based on Fig. S4, 15% appears to be optimal/consistent with the experiment, rather than 10%. It would be useful to\ncomment on this difference.\nAuthors’ response: We thank the reviewer for raising this important point. We assessed the convergence of equilibrium\nproperties for the disordered protein fragments by systematically analyzing their secondary structure propensities at different\nintervals throughout the simulations.\nAs detailed in the Methods section and explicitly shown in Figure S11, the residual helix fraction calculated from the entire\ntrajectory was consistent with values obtained after excluding the initial 50 ns and 100 ns of simulation data. This\nconsistency demonstrates that our PT-WTE simulations achieved converged equilibrium sampling of secondary structure.\nTo clearly communicate this evidence to readers, we have revised the manuscript sentence as follows:\nRevised sentence:\n“...For these fragments, we employed PT-WTE simulations to ensure converged sampling of their equilibrium properties,\nas demonstrated by our secondary structure convergence analysis (see Methods and Figure S11).”\nAdditional details on scaling factor testing are provided in our response to the related minor comment below.\n(5) Fig. 4 — It is not clear where the experimental data shown in this figure comes from. Ref. 59 cited in the paragraph\nabove Fig. 4 does not show this data.\nAuthors’ response: We thank the reviewer for pointing out this ambiguity. The experimental data in Figure 4 are derived\nfrom temperature-dependent carbonyl chemical shift measurements of (AAQAA)₃ reported by Shalongo and Stellwagen\n(Ref. 59). The helical fractions were calculated using the chemical shift data in their Table 1 and Equation 2, which relates\nchemical shift deviations to helical population.\nWe have added clarification to the figure caption as follows:\n“Figure 4. Temperature dependence of structural properties of the (AAQAA) model peptide. Helix fraction of (AAQAA)\n3 3\nas a function of temperature simulated in ff03ws, ff03w-sc, and ff99SBws, respectively. Experimental data (Exp.) were\ncalculated from the carbonyl chemical shifts reported by Shalongo and Stellwagen.59 The inset shows the residual helix\nfraction at 300 K. Error bars for the simulation data indicate the standard deviation.”\n(6) line 297-298 — “We strategically selected 44-residue fragments of these proteins to optimize computational efficiency\nwhile maintaining biological relevance and employed PT-WTE simulations to achieve converged equilibrium properties.”\nWhere is the convergence of the equilibrium properties shown?\nAuthors’ response: We thank the reviewer for raising this important point. We assessed the convergence of equilibrium\nproperties for the disordered protein fragments by systematically analyzing their secondary structure propensities at different\nintervals throughout the simulations.\nAs detailed in the Methods section and explicitly shown in Figure S11, the residual helix fraction calculated from the entire\ntrajectory was consistent with values obtained after excluding the initial 50 ns and 100 ns of simulation data. This\nconsistency demonstrates that our PT-WTE simulations achieved converged equilibrium sampling of secondary structure.\nTo clearly communicate this evidence to readers, we have revised the manuscript sentence as follows:\nRevised sentence:\n“...For these fragments, we employed PT-WTE simulations to ensure converged sampling of their equilibrium properties,\nas demonstrated by our secondary structure convergence analysis (see Methods and Figure S11).”\nAdditionally, please refer to our response to the related minor comment below regarding the rationale for selecting the\nspecific 44-residue fragments.\n(7) lines 302-303 — “Figure 6A shows excellent agreement between simulated and experimental data in the helix-forming\nregion (residue 321–330) in TDP-43310–350” -> this statement implies that the simulation results agree with experiment\nfor all data points, which does not match what is shown in the figure.\nAuthors’ response: We thank the reviewer for this comment, which helps improve the precision of our language. We agree\nthat “excellent agreement” is an overstatement given the quantitative differences observed for some residues. Our intention\nwas to emphasize that both force fields captured the experimentally observed helix propensity specifically within the known\nhelix-forming region of TDP-43 (residues 321–330), rather than suggesting agreement across the entire sequence.\n310-350\nWe have revised the sentence to reflect the correct prediction of helical propensity of this fragment.\n“Figure 6A shows that both force fields successfully captured the experimentally observed helix propensity within the core\nhelix-forming region (residues 321–330) in TDP-43 , while deviations are present in other regions of the sequence.”\n310–350\n(8) lines 306-307 — “The RMSDs of simulated chemical shifts relative to experimental measurements for both force\nfields were less than the prediction error of SPARTA+ (± 1 ppm for carbon chemical shifts66), highlighting their\nreasonable accuracy in modelling IDP conformational ensembles.” - It would be useful to provide the data to support this\nstatement.\nAuthors’ response: We thank the reviewer for the suggestion. The RMSD values referenced in this sentence were computed\nas part of our secondary chemical shift analysis and are shown in Figures 6 and S6. For each IDP fragment, we calculated\nthe RMSDs of the simulated ΔδC − ΔδC values relative to the experiment. These RMSDs were consistently below 1 ppm,\nα β\nwhich is within the reported prediction error of SPARTA+ for carbon chemical shifts (±1 ppm). To improve clarity, we\nhave revised the sentence to explicitly state this and added a reference to the supporting figures.\n“The RMSE values of simulated chemical shifts (ΔδC − ΔδC ) relative to experimental values for both force fields were\nα β\nconsistently below 1 ppm (Figures 6 and S6), which is within the prediction error of SPARTA+,70 implying the reasonable\naccuracy of both force fields in modeling IDP structural ensembles.”\n(9) lines 316-317 — “In addition, we compared secondary structure propensities calculated from simulation trajectories\nusing DSSP 67” —> DSSP as described in ref. 67 does not provide secondary structure propensities. Further details on how\nthe authors did this are needed. It would also be useful to explain why the authors chose to compare helix fraction predicted\nfrom NMR chemical shifts using the delta2d web server rather than directly comparing to the experimental chemical shifts.\nAuthors’ response: We thank the reviewer for this thoughtful comment. The reviewer is correct that the original DSSP\nalgorithm (as described in Ref. 67) provides per-residue, per-frame secondary structure assignments rather than directly\ncalculating propensities. In our analysis, we used the GROMACS tool “gmx do_dssp” to assign DSSP-defined secondary\nstructures at each trajectory frame, then calculated helix fractions by averaging the number of frames in which each residue\nwas assigned to a helical state. To clarify this, we have updated the Methods section to explicitly describe how secondary\nstructure propensities were derived from DSSP assignments.\n“...Secondary structure assignments from simulations were determined using “gmx do_dssp”,95 which implements the DSSP\nalgorithm.67 The propensity for a given structure (e.g., helix fraction) for each residue was determined by averaging the\nper-frame secondary structure assignments over the entire trajectory…”\nRegarding the use of the δ2D web server, we chose this approach because it enables direct inference of per-residue secondary\nstructure populations (helix, sheet, and coil) from experimental NMR chemical shifts. While direct comparison to chemical\nshift values (e.g., C , C ) is informative, δ2D provides a more interpretable and structure-specific comparison for\nα β\nbenchmarking simulation-derived secondary structure distributions. We have added a sentence to the manuscript to clarify\nthis choice.\n“In addition, we calculated secondary structure propensities from the simulation trajectories based on the DSSP algorithm71\nfor each frame and averaged the helical assignments per residue over time. These propensities were compared to those\npredicted from experimental NMR chemical shifts using the δ2d web server,72 which converts chemical shifts into residue-\nlevel secondary structure populations (Figure 6C, D).”\n(10) The authors present several new results in Discussion/Conclusions that would be more appropriate in the Results\nsection (Fig. S9-S10).\nAuthors’ response: We thank the reviewer for this helpful structural suggestion. We agree that the analyses presented in\nFigures S9 and S10 should be included in the Results section rather than the Discussion. We have revised the manuscript\naccordingly by moving these analyses and their interpretations to improve the logical flow:\n- The NMR-derived order parameter analysis for ubiquitin (Figure S9) has been moved to the end of the “Modifications to\nthe protein-water scaling scheme in amber ff03ws” section.\n“To further assess force field performance on folded proteins, we investigated their ability to capture backbone dynamics\non the ps-ns timescale by calculating the amide order parameters (S²) for ubiquitin and compared them to the NMR-derived\nvalues (Figure S7). While the refined ff03w-sc force field successfully stabilized the overall structure of ubiquitin, it\nexhibited the known tendency of ff03-based force fields to overestimate backbone flexibility, particularly in loop regions\nwhere simulated S² values were consistently lower than experimental measurements (Figure S7A). In contrast, ff99SB-\nbased force fields demonstrated better performance, yielding S² order parameters in excellent agreement with NMR across the\nentire protein and indicating more accurate descriptions of both rigid and flexible regions (Figure S7B). These results\nhighlight the differential performance of the two force field families in capturing local backbone dynamics.”\n- The residue-type–averaged chemical shift deviations for low-complexity fragments (Figure S10) are now described in a\nnew paragraph at the end of the “Secondary structure predictions for IDPs” section.\n“To further evaluate the force fields for systematic, residue-specific biases, we analyzed the chemical shift deviations\naveraged by amino acid type across all four low-complexity fragments (Figure S14). This analysis shows that for most\nresidue types, both force fields yield average deviations close to zero, indicating good overall accuracy. However, some\ndifferences emerge for polar residues. The ff03w-sc force field shows notable positive deviations for serine and threonine\n(similar to ff03ws64), suggesting a persistent bias toward helical conformations for these residue types. These deviations\nare substantially reduced in the ff99SBws-STQ force field, which demonstrates the success of its targeted torsional\nrefinements for serine, threonine, and glutamine and highlights its improved accuracy for modeling polar-rich disordered\nsequences.”\nThe figures in the Supplementary Information have been reordered accordingly; the former Figures S9 and S10 now appear\nas Figure S7 and Figure S14, respectively.\nThese revisions ensure that all primary results are presented in the Results section, with the Discussion now focused on\ninterpretation and broader implications.\n(11) Line 482 — The link to the force field files does not work.\nAuthors’ response: We thank the reviewer for bringing this to our attention. We have fixed the broken link, and the force\nfield files are now accessible at the updated URL:\nhttps://bitbucket.org/jeetain/all-atom_ff_refinements.\nWe have also verified that the link works correctly and updated the manuscript accordingly.\nMinor comments:\n- To enhance the clarity of the results, the authors should state explicitly which scaling factor for the side chains is being\nused at different points in the manuscript. This becomes confusing because different possible scaling factors have been\nconsidered.\nAuthors’ response: We thank the reviewer for this suggestion. The 15% scaling factor shown in Figure S4 was tested solely\nto demonstrate the trade-off between IDP conformational dimension and folded protein stability, but 10% scaling was\nselected as the final choice for our final force field.\nTo avoid confusion, we have revised the relevant paragraph in the Results section to clearly distinguish between exploratory\ntests (e.g., the 15% scaling) and the definitive 10% scaling used in the final ff03w-sc force field throughout all main\nsimulations and analyses.\n“We note that a 10% increase in protein side chain-water interactions used in ff03w-sc represents a near-optimal balance\nbetween maintaining folded protein structural stability and preventing over-collapse of disordered proteins. To further\nillustrate this trade-off, we also tested a stronger 15% scaling factor. While this higher affinity can yield more realistic\nchain dimensions for disordered proteins, it also leads to destabilization of folded proteins (Figure S8A-C). Based on this\nanalysis, we adopted the 10% scaling factor as the final choice for ff03w-sc throughout this study. Notably, applying a side\nchain scaling of 10% to ff99SBws resulted in overly collapsed conformations of CspM34 (Figure S8D), highlighting the\nsensitivity of this force field to scaling adjustment.”\nPlease also see our response to the specific comment (4) above regarding the choice of the 10% scaling factor.\n- Line 257 — Should be Figure 4, not Figure 3.\nAuthors’ response: We thank the reviewer for catching this error. We have corrected the figure reference from Figure 3 to\nFigure 4 in line 257.\n- Line 274-278 — Reference should be provided for this sentence.\nAuthors’ response: We thank the reviewer for pointing out this missing citation. We have added the appropriate reference\nfor this sentence in the revised manuscript.\n- Line 296-298 — How were the stretches of 44 residues chosen?\nAuthors’ response: We thank the reviewer for this question regarding our choice of fragments. As stated in the manuscript,\nthe selected lengths (44 residues) for TDP-43, FUS, and RNA Pol II fragments were chosen to balance computational\nefficiency and biological relevance. Specifically, these segments encompass key regions that are well-characterized\nexperimentally for their conformational propensities and secondary structure elements. By selecting these representative\nstretches, we ensured adequate conformational sampling and reliable comparison with experimental data, enabling robust\nvalidation of our force fields while maintaining manageable computational demands. We have clarified this rationale in the\nmanuscript accordingly.\n“We strategically selected 44-residue fragments of these proteins to optimize computational efficiency while focusing on\nbiologically relevant regions for which detailed NMR chemical shifts and secondary structure data are readily available,\nallowing direct validation of our simulation results. For these fragments, we employed PT-WTE simulations to ensure\nconverged sampling of their equilibrium properties, as demonstrated by our secondary convergence analysis (see Methods\nand Figure S10).”\n- Line 322 — Should be Figure 6c,d not 5c,d\nAuthors’ response: We thank the reviewer for careful proofreading. We have corrected the figure reference from Figures 5c\nand 5d to Figures 6c and 6d in line 322.\n- Fig S3 — “ff99SBwc” in y-axis labels should be “ff99SBws”.\nAuthors’ response: We thank the reviewer for identifying this typo. We have corrected the y-axis label from “ff99SBwc”\nto “ff99SBws” in the revised Figure S3.\nRESPONSE TO REVIEWERS’ COMMENTS\nReviewer #2 (Remarks to the Author):\nThe authors have addressed the majority of my earlier comments in their revised manuscript. There are a few\nremaining points to address that could clarify and strengthen their manuscript:\n(1). Lines 237-238\n“Based on this analysis, we adopted a 10% scaling as the final choice for ff03w-sc throughout this study.”\nThe authors have partially explained why they adopted this choice of scaling. However, it should be noted that\nthis choice does not consistently result in ensembles that are in closest agreement with the experimental data (e.g.,\nin Figs. 3A, S8 A&D).\nAuthors’ response: We agree that a 10% sidechain-water scaling does not always produce the per-system best fit\nto every observable. Our aim was transferability, an optimal scaling that preserves folded-state stability while\nappropriately expanding disordered ensembles across systems, rather than optimizing any one metric in isolation.\nIn our tests, stronger scaling sometimes improved IDP dimensions but compromised folded-protein stability,\nwhereas 10% provided a robust cross-system balance. We have clarified this rationale in the manuscript at Lines\n237–239.\n“We note that 10% increase in protein side chain-water interactions used in ff03w-sc represents a near-optimal\nbalance between maintaining folded protein structural stability and preventing over-collapse of disordered\nproteins. To further illustrate this trade-off, we also tested a stronger 15% scaling factor. While this higher affinity\ncan yield more realistic chain dimensions for disordered proteins, it also leads to destabilization of folded proteins\n(Supplementary Figures 8A-C). Therefore, although 10% scaling may not always yield the closest agreement\nwith experimental data for every individual system, it provides robust transferability across diverse protein\nsystems. Based on this analysis, we adopted 10% scaling as the final choice for ff03w-sc throughout this study.”\n(2) Lines 325-329\n“Additionally, 3JHN-HA scalar couplings calculated for the RS-peptide showed improved agreement with\nexperiment relative to ff03w-sc and ff99SBws (Figure S9F, Table S5), implying an improvement in local\nbackbone configurational sampling. These results demonstrate the capability of ff99SBws-STQ′ to capture both\nglobal dimensions and local structural properties across a wide range of protein systems.”\nIt is unclear what the authors mean by improved agreement here, especially for ff03w-sc. The results shown in\nFig. S9F (left) are not in agreement with the experimental data, even qualitatively. The authors acknowledge this\ninconsistency in their response to the reviewers, and the details provided in the response would be beneficial to\nadd to the manuscript.\nThe details I am referring to are:\n“Across all force fields, we evaluated couplings using multiple Karplus parameter sets, including DFT1,\nDFT2, and parameters from Vogeli et al. and Larson et al., as well as the original Karplus form. The experimental values\nare presented as a gray-shaded band (mean ± SD), facilitating direct visual comparison. As seen in the figure:\n• ff99SBws-STQ′ shows the closest agreement with experiment across the entire RS-rich sequence, with values\nfalling well within experimental uncertainty for most residues, particularly when using the DFT1 or DFT2 Karplus\nparameters.\n• ff99SBws also performs reasonably well, though it slightly underestimates the couplings in some central serine\nresidues.\n• ff03w-sc only captures the qualitative trend while exhibiting systematic underestimation of couplings, with\nlarger deviations across most positions.\nWe note that this behavior of ff03w-sc is consistent with previously published results for ff03w and ff03ws,\nboth of which show similar underestimation of scalar couplings for RS-rich sequences (see Ref. 18, Figure 5A). As\nnoted in Ref. 18, this suggests that the underlying backbone φ dihedral angle preferences of the ff03-based force\nfield family deviate from experiment for this peptide, a behavior that is largely independent of the specific water-\nprotein scaling scheme. Given that ff03w-sc inherits its backbone parameters from ff03w, the observed behavior\nis in line with previous work. Overall, these results support the conclusion that ff99SBws-STQ′ provides improved\nlocal conformational sampling for RS-rich sequences, and they further validate the force field refinement against\nindependent NMR observables.”\n“These results demonstrate the capability of ff99SBws-STQ′ to capture both global dimensions and local\nstructural properties across a wide range of protein systems.”\nThis is an overstatement (specifically, local structural properties) given that scalar couplings were compared for\nonly the RS peptide and Ala5 (and the results were not consistent with the experimental data).\nAuthors’ response: We appreciate the opportunity to improve our manuscript’s clarity. For the RS peptide, we\nevaluated 3J couplings using multiple Karplus parameterizations (original Karplus, Vogeli et al., DFT1,\nHN-HA\nDFT2, and Larson et al.) and compared them to experiment, shown as a gray band (mean ± SD) in Supplementary\nFigure 9F (with numerical values in Supplementary Table 5). Under these conditions, ff99SBws-STQ′ yields the\nclosest overall agreement, with most residues falling within experimental uncertainty—particularly with\nDFT1/DFT2. ff99SBws performs reasonably but tends to underestimate couplings for the central serines. By\ncontrast, ff03w-sc exhibits a systematic underestimation across much of the sequence and only partial qualitative\nagreement, consistent with prior reports for ff03-family models on RS-rich sequences (Rauscher et al., 2015). To\navoid overstating the scope, we have revised the text to (i) specify the Karplus sets and how agreement was\nassessed, and (ii) temper the concluding claim about local structure to the systems examined here, noting that\nbroader generality will require additional test cases.\n“Additionally, we evaluated 3J scalar couplings for the RS-peptide using multiple Karplus parameter sets\nHN-HA\n(original Karplus,70 Vogeli et al.,71 DFT1,72 DFT2,72 and Larson et al.;73 Supplementary Figure 9F and\nSupplementary Table 5). ff99SBws-STQ′ shows the closest overall agreement with experimental values, with most\nresidues falling within experimental uncertainty (especially with DFT1/DFT2), whereas ff99SBws slightly\nunderestimates couplings for some central serine residues. In contrast, ff03w-sc systematically underestimates\ncouplings across most sequence positions, consistent with prior observations for ff03-based models.18 Together\nwith global-dimension benchmarks, these results demonstrate that ff99SBws-STQ′ accurately captures global\ndimensions across tested systems and improves local backbone sampling for the RS-peptide. However, a broader\nassessment of local structure across additional sequences remains an important direction for future work.”\n(4) Lines 393-394\n“Collectively, these Ala5 simulations provide additional quantitative validation for both force fields.”\nSimilar to the comment above, it is an overstatement to refer to this as “validation,” given the inconsistency with\nthe experimental data.\nAuthors’ response: We thank the reviewer for this suggestion. We agree that “validation” overstates the\nconclusions given the modest deviations from experimental data, particularly for ff99SBws-STQ'. Our intent was\nto use Ala as a quantitative assessment of local backbone sampling rather than a definitive validation. We have\n5\nrevised the manuscript to use more appropriate terminology.\n“Collectively, these Ala simulations provide an additional quantitative assessment for both force fields. Despite\n5\ntheir differing refinement strategies, both force fields exhibit reasonable agreement with experimental scalar\ncoupling data, thereby supporting their ability to achieve a reliable sampling of the conformational landscape for\nsmall peptides.”\n(5) Lines 446-447\n“Both force fields preserved high fractions of native interfacial contacts (>95%) across all complexes, consistently\nmaintaining key interfacial interactions.”\nThis statement is inconsistent with Fig. S17. It implies that the fraction of interfacial native contacts is always\ngreater than 95%.\nAuthors’ response: We thank the reviewer for their careful attention to detail. We agree that our original statement\nimplied a strict lower bound (>95%), which was not consistently supported by the time-evolution traces shown\nin Supplementary Figure 17. We have revised the sentence to more accurately reflect the data. The new phrasing\nremoves the strict “>95%” threshold and instead describes the performance in terms of the median values\nobserved across the simulations, which still conveys the excellent stability of the interfaces without\nmisrepresenting the data.\n“Both force fields preserved a high fraction of native interfacial contacts, with median stability values near 95%\nacross the simulations.”\n(6) Lines 411-413\n“As shown in Figures 7 and S15, ff03w-sc reasonably reproduces the experimental melting behavior for both\nmini-proteins, which are consistent with previous observations reported for ff03ws. Folded fractions closely\ntrack the experimental curves across the temperature range.”\nThis description is an overstatement given the clear inconsistency between the experiment (black line) and\nsimulation (red points) in Fig. 7. Lines 539-544 in the discussion would also need to be modified accordingly,\nsuch that the findings are not overstated.\nAuthors’ response: We thank the reviewer for this suggestion. Our intent was to convey that ff03w-sc yields a\nhigher folded population and captures the qualitative melting trend for Trp-cage and chignolin, not that it\nreproduces the quantitative experimental curves across the full temperature range. We have revised the language\naccordingly and now explicitly note the quantitative deviations.\nLines 419-421:\n“As shown in Figure 7 and Supplementary Figure 15, ff03w-sc captures the qualitative cooperative melting\nbehavior for both mini-proteins, consistent with previous observations reported for ff03ws,25 although\nquantitative deviations from experiment are evident across the temperature range. In contrast, ff99SBws-STQ′\nsignificantly underestimates the folded populations of both systems, particularly near physiological temperatures,\nindicating destabilization of native conformations.”\nLines 547-549:\n“Furthermore, folding simulations of mini-proteins revealed that ff03w-sc yields higher folded populations and\ncaptures the qualitative melting transition for both Trp-cage and chignolin, despite quantitative deviations from\nexperiment, while ff99SBws-STQ′ significantly underestimated folded fractions across the temperature range, a\nbehavior also observed with other force fields such as ff99SB-disp and charmm36m.29 This suggests that\naccurately balancing local backbone propensities with global folding cooperativity remains a central challenge\nfor developing widely applicable force fields.”\nAdditionally, we have revised the conclusion to better acknowledge the inherent limitations of force field\ndevelopment and emphasize the need for continued methodological improvements.\n“In conclusion, this study advances the ongoing efforts to refine protein force fields, enhancing their utility in\nbiomolecular simulations. The proposed force fields— ff03w-sc and ff99SBws-STQ—achieve a balanced\ndescription of both folded and IDP ensembles and are well-suited for studying a diverse range of systems,\nincluding folded proteins, disordered regions, and multiprotein complexes. These advancements help narrow the\ngap between theoretical predictions and experimental observations, enabling more predictive simulations of\nprotein behavior. However, it is important to recognize that no single force field can perfectly capture the\ncomplexity of all protein systems—each model inevitably has distinct strengths and limitations depending on the\nspecific structural and dynamic properties being investigated. Given this fundamental challenge, future work\nshould focus on developing strategies to further refine these force fields to better balance local structural\npropensities with global folding cooperativity. Understanding the inherent trade-offs in force field design will be\ncrucial for guiding these improvements and determining optimal applications for each model. Such efforts will\npave the way for extending these improvements to larger, multidomain proteins and exploring their applicability\nunder biologically relevant conditions, such as crowded intracellular environments and systems with post-\ntranslational modifications.”"}]

Decision letter