Structure-guided forcefield optimization

Accurate modeling of biomolecular systems requires accurate forcefields. Widely used molecular mechanics (MM) forcefields obtain parameters from experimental data and quantum chemistry calculations on small molecules but do not have a clear way to take advantage of the information in high-resolution macromolecular structures. In contrast, knowledge-based methods largely ignore the physical chemistry of interatomic interactions, and instead derive parameters almost exclusively from macromolecular structures. This can involve considerable double counting of the same physical interactions. Here, we describe a method for forcefield improvement that combines the strengths of the two approaches. We use this method to improve the Rosetta all-atom forcefield, in which the total energy is expressed as the sum of terms representing different physical interactions as in MM forcefields and the parameters are tuned to reproduce the properties of macromolecular structures. To resolve inaccuracies resulting from possible double counting of interactions, we compare distribution functions from low-energy modeled structures to those from crystal structures. The structural and physical bases of the deviations between the modeled and reference structures are identified and used to guide forcefield improvements. We describe improvements resolving double counting between backbone hydrogen bond interactions and Lennard-Jones interactions in helices; between sidechain-backbone hydrogen bonds and the backbone torsion potential; and between the sidechain torsion potential and Lennard-Jones interactions. Discrepancies between computed and observed distributions are also used to guide the incorporation of an explicit Cα-hydrogen bond in β sheets. The method can be used generally to integrate different sources of information for forcefield improvement.


INTRODUCTION
High-resolution protein structures provide an invaluable source of information for forcefield testing and improvement. The parameters of knowledge-based energy functions are entirely derived from protein structures. With the assumption that distributions of interatomic distances and other protein structure properties are independent and Boltzmann distributed, the underlying energy is obtained simply by computing the negative logarithm of the observed distributions. 1 However, the distributions are far from independent, and this approach also has the disadvantage of neglecting the large body of knowledge on the physical chemistry of interatomic interactions. 2 Molecular mechanics (MM) forcefields on the other hand derive most parameters from experimental and quantum chemistry data on small molecule data rather than proteins. 3-8 The total system energy is expressed as the sum of terms with clear physical origins. The applicability and transferability of such forcefields to macromolecules has been demonstrated by showing that native structures are at least metastable in molecular dynamics (MD) simulations and by reproducing vibration spectroscopic data. 6,8 However, these tests only probe the accuracy of a forcefield in the neighborhood of the native state; longer time scale simulations are now starting to test forcefield accuracy over a larger range of conformational space. 9,10 The Rosetta forcefield, like MM forcefields, expresses the total system energy as the sum of physically interpretable terms but uses protein structure data both to refine parameters and to model contributions that are difficult to obtain by other means. 11 For example, the Rosetta forcefield supplements Lennard-Jones and implicit solvation terms used in MM forcefields with protein structure-derived sidechain and backbone torsion potentials.
Double counting of the same physical interaction by two different forcefield terms can result in overall forcefield inaccuracies. It is not straightforward to systematically correct this problem. In MM forcefields, the backbone torsion potential, combined with the rest of the forcefield, may incorrectly bias the balance between helix and sheet structures. This problem has been addressed by compar-Additional Supporting Information may be found in the online version of this article.
ABSTRACT Accurate modeling of biomolecular systems requires accurate forcefields. Widely used molecular mechanics (MM) forcefields obtain parameters from experimental data and quantum chemistry calculations on small molecules but do not have a clear way to take advantage of the information in high-resolution macromolecular structures. In contrast, knowledge-based methods largely ignore the physical chemistry of interatomic interactions, and instead derive parameters almost exclusively from macromolecular structures. This can involve considerable double counting of the same physical interactions. Here, we describe a method for forcefield improvement that combines the strengths of the two approaches. We use this method to improve the Rosetta all-atom forcefield, in which the total energy is expressed as the sum of terms representing different physical interactions as in MM forcefields and the parameters are tuned to reproduce the properties of macromolecular structures. To resolve inaccuracies resulting from possible double counting of interactions, we compare distribution functions from low-energy modeled structures to those from crystal structures. The structural and physical bases of the deviations between the modeled and reference structures are identified and used to guide forcefield improvements. We describe improvements resolving double counting between backbone hydrogen bond interactions and Lennard-Jones interactions in helices; between sidechain-backbone hydrogen bonds and the backbone torsion potential; and between the sidechain torsion potential and Lennard-Jones interactions. Discrepancies between computed and observed distributions are also used to guide the incorporation of an explicit Ca-hydrogen bond in b sheets. The method can be used generally to integrate different sources of information for forcefield improvement.
ing molecular dynamics simulations to NMR coupling and relaxation data and quantum chemistry calculations 12,13 and adjusting the torsional potentials accordingly. For residue-resdiue knowledge-based potentials, an iterative improvement approach has been described 14 that compares computed and experimentally observed minima.
In this study, we use an iterative approach to detect and remedy problems resulting from double counting in forcefields based on the properties of energy minima distributed throughout conformational space. Systematic structural comparisons between the X-ray crystal structures and refined Rosetta models generated both near and far from their native state are used to track down errors in the standard Rosetta forcefield. The physical origin of the errors are identified and used to guide correction of the individual forcefield terms and to explicitly remove double counting.

Crystal structure dataset
A set of high-resolution crystal structures was used as the reference for atom-atom distance, angle, and torsion distribution calculations. X-ray crystal structures of 1257 proteins were selected using PISCES 15,16 with resolution better than 1.5 Å , crystallographic R factor better than 0.3, and maximum sequence identity of 25%. Water and ligands were removed.

Energy landscape generation
Energy landscape calculations and characterization are described in detail by Tyka et al. 17 In brief, low-resolution models spanning a broad range of RMSD to the native structure (0$20 Å ) are generated using the Rosetta ab initio folding protocol 18 for 110 proteins. These 110 proteins are selected to include a variety of secondary structural elements and structural features. They include 24 all-alpha, 29 all-beta, and 57 alpha-beta proteins. Among these proteins, 17 of the structures bind a ligand, 60 are oligomeric, and 37 contain disulfide bonds.
The Rosetta full atom refinement protocol 18,19 is then applied to search for local minima in the vicinity of each low-resolution model with either the standard or the optimized energy function. For each protein, 100,000 allatom refined models are generated, and the models are then placed into bins based on their RMSD to the native structure (bin width 0.5 Å ). In each RMSD bin, the 20% of structures with lowest energies are collected for distribution calculations. This ensures that the distributions reflect contributions spanning a large range of structures not just conformations near the native structure.
The energy gap between native and non-native structures is taken to be the (average energy of the lowest 1% of structures that are less than 2 Å RMSD from the native structure) minus the (average energy of the lowest 1% of structures that are greater than 4 Å from the native structure).

Distribution calculation
Atom-atom radial, angular, and dihedral distribution functions were collected for the 1257 protein crystal structure dataset and the low-energy computed structures for the 110 selected proteins described in the previous section. The low-energy Rosetta models from different RMSD bins were pooled together in the distribution calculation. Backbone atom-atom radial distributions were determined between all backbone atom pairs for each secondary structure type (secondary structure was designated as a-helix, b-strand, or loop as determined by the DSSP algorithm 20 ). From crystal structures, there are a total of 2.7 3 10 5 helix, 0.6 3 10 5 b-sheet, and 1.0 3 10 5 loop residues. Atom pairs less than 10 Å apart were collected for 1.2 3 10 6 helix-helix residue pairs, 0.7 3 10 6 strand-strand pairs, and 0.9 3 10 6 loop-loop pairs. In each iteration of Rosetta modeling, distributions were computed from 1.3 3 10 7 helix, 1.0 3 10 7 b-sheet, and 1.3 3 10 7 loop residues. The atom-atom radial distribution function is the average density of an atom2 at a distance r from atom1, qðr; atom2jatom1; ss1; ss2Þ ¼ N ðr; atom1; atom2; ss1; ss2Þ 4pr 2 Dr3N ðatom1; ss1; ss2Þ ð1Þ where the bin width Dr 5 0.1 Å bin, N ðr; atom1; atom2jss1; ss2Þ is the smoothed [Eq. (2)] counts of atom pairs (atom1 and atom2) at distance r with secondary structure ss1 and ss2, and N ðatom1jss1; ss2Þ is the sum of N(r,atom1, atom2|ss1, ss2) over all atom2 types.
The counts N ðr; atom1; atom2jss1; ss2Þ were smoothed using a Gaussian kernel: N ðr i ; atom1; atom2; ss1; ss2Þ where N o ðr j ; atom1; atom2; ss1; ss2Þ is the total counts at r j before smoothing, and the bin width and variance s are 0.1 Å . Distributions were also determined for angles between hydrogen-bonding atoms. Two angles are measured: Y, the angle formed by the donor heavy atom-donor pro-ton-acceptor atom triplet, and X the angle formed by the donor proton-acceptor atom-acceptor base atom triplet. 21 A total of 5.7 3 10 4 helix and 3.7 3 10 4 b-sheet backbone hydrogen bonds were collected from crystal structures, and 5 3 10 6 were collected from low-energy Rosetta models in each iteration of forcefield optimization. Angular distributions were calculated using: where Dy, the bin width, is 58, N y, total is the total number of hydrogen bonds, N y, i is total counts in the angle bin y i . Counts were first smoothed using Gaussian kernel smoothing with variance 58. Backbone dihedral distributions (Ramachandran distributions) were collected for each residue type in each secondary structure. In the reference crystal structure set, for a given residue type and secondary structure, the numbers of residues collected for the distribution calculations range between 2,000 and 5,000. In the low-energy Rosetta models, the numbers of residues collected for each residue type and secondary structure combination range between 2 3 10 5 and 4 3 10 5 . Counts were binned into two-dimensional 108 bins in u and w and then smoothed using a kernel width of 108.
For sidechain distributions, we followed previous work 22-24 and classified sidechain conformations into discrete rotamer bins based on their v angles and tabulated statistics separately by 108 backbone u/w bins. Only rotamers with less than 20% accessible surface area are counted in the distribution. The rotamer definitions and the treatment of the terminal v angle for Asn, Asp, Gln, Glu, His, Phe, and Tyr (5 or 10 degree bins) were as in the 2008 Dunbrack rotamer library. 24 Gaussian kernel smoothing was applied to both the backbone dependency and the terminal v angles for the sidechains listed above; for the remaining fully rotameric sidechains, no smoothing was applied.

Identification of inaccurately modeled features
We assess the forcefield at each iteration by computing the difference in the logarithms of the distributions of each of the above features between the reference crystal structures and low-energy Rosetta structures.
where q i,sample and q i,ref are the smoothed Rosetta model and crystal structure distributions. Large differences between the crystal structure distribution and modeled structures indicate potential areas for forcefield improvement. As we wish to focus improvements in frequently sampled regions and to reduce noise in the log-difference due to small counts, we scale the difference by the average frequency with which each bin is populated: where q mean is the Boltzmann average probability over all bins in the crystal structures: We focused subsequent analysis on regions where the scaled log difference between the distributions was greater than 1.

Forcefield optimization
The regions where large differences are found between Rosetta modeled structures and crystal structures were analyzed, and the physical origin of the discrepancies were identified to determine the potential term that needs improvement. Identification involved inspection of the structural contexts of the discrepancies and analysis of the potential terms contributing to the Rosetta structure distributions.
Once the potential term requiring correction was identified, an iterative approach was applied to update the form of the potential guided by the difference between Rosetta models and crystal structures. For the Ca hydrogen bond potential (described below), Ramachandran potential and rotamer potential, the collected distribution function has the same dimension as the potential function, and therefore, the log difference between the Rosetta model and crystal structure distributions [Eq. (5)] can be directly subtracted from the potential term at each iteration: where E i (n) is the potential function used in the nth iteration, and De i (n) is the log difference in Eq. (5) at the nth iteration. The new potential is then used for a new round of full atom refinement to generate a new set of near-native and non-native structures. Iterations of structure model generation, distribution calculation, and potential correction continue until De i (n) is below 0.5 for all bins. It took four iterations for the Ca hydrogen bond potential to reach convergence, 12 iterations for Ramachandran potential, and eight iterations for rotamer potential. A similar iterative approach was used previously to improve a simple pairwise residue-residue potential. 14 The hydrogen bond potential is the sum of energy terms based on hydrogen bond distance and angles, where r is the distance between the donor proton and acceptor atom, y is the hydrogen bond angle formed by donor heavy atom-donor proton-acceptor atom, and x is formed by donor proton-acceptor atom-acceptor base atom as described earlier. 21 At each iteration, the potential forms of e Y and e X are maintained, but the peak of the potential is shifted: where DY n is the difference in the peak position between the input e Y (y) potential and the sampled distribution in the nth iteration. A similar correction is applied to e X (x). The update of a single parameter rather than the entire potential function reduced over-fitting artifacts.
The new Ca hydrogen bond potential was defined simply as a distance-dependent interaction between Ha and O. The starting guess at the potential was simply the log of the crystal structure HaÀ ÀO distance distribution function, and then the iterative approach described above was applied to optimize the potential until the modeled HaÀ ÀO radial distribution function matched that of crystal structures.

Independent benchmark test
An additional benchmark test was applied to a dataset independent from the 110 proteins used for the optimization. This benchmark test uses 55 protein structures from the CASP8 experiment. 25 The HHSearch protocol 26 was used to generate alignments to the pre-CASP8 database of template structures. Complete models based on those alignments were generated using Rosetta loop modeling. 27,28 Then, the Rosetta full atom refinement protocol 18,19 was applied to search for low-energy models with either standard or the optimized energy function. For each protein, 10,000 all-atom refined models were generated.

RESULTS
As described in the Methods section, we generated large ensembles of conformations sampling local minima throughout the energy landscape for a set of 110 proteins of known structure. From these ensembles, we computed distributions of interatomic distances and bond torsion angles. These distributions were compared with those observed in high-resolution protein structures, and the atom pairs and torsions for which the distributions in the computed energy minima differed significantly with crystal structures were flagged. For the majority of atom pairs and torsion angles, the distributions match quite closely; several examples are shown in Figures 1 and 2.
Significant discrepancies between the distributions in crystal structures and computed energy minima likely reflect errors in the forcefield. Inspection of these discrepancies suggested that the majority of the errors reflect double counting of the same physical chemistry by two different energy function components. Other discrepancies were traced to missing energetic contributions. These errors and their resolution are described in the following sections.

Helix backbone hydrogen bonding
The distribution of the hydrogen bond distance between the backbone nitrogen and the carbonyl oxygen is similar in crystal structures and low-energy Rosetta models (Fig. 1). However, the distribution of distances between pairs of carbonyl oxygen atoms differs between the two sets of structures. As shown in Figure 1(A), the distribution peak at 3.4 Å in crystal structures is shifted to 3.6 Å in Rosetta models. This peak is due to pairs of backbone oxygens in residues adjacent in the sequence [ Fig. 1(E)]. These atoms primarily interact with each other in the Rosetta forcefield via Lennard-Jones interactions. The helix hydrogen bond angular distributions for the angles formed by the atom triplets (NÀ ÀH)À ÀO and HÀ À(O¼ ¼C) also differ between Rosetta energy minima and crystal structures. The peak for the (NÀ ÀH)À ÀO angle is shifted from 1628 in crystal structures to 1578 in Rosetta models, and the HÀ À(O¼ ¼C) angle is shifted from 1528 to 1458.
To identify which energy terms contribute to the observed shifts of the peaks in the hydrogen bond angle distributions, low-energy models were generated with all Rosetta energy terms up-or down-weighted individually. When the Lennard-Jones repulsion term was changed, changes in hydrogen bond angle peaks were the most significant, suggesting the shift of the peak is likely due to the coupling between the Lennard-Jones term and hydrogen bond energy: as neighboring oxygens are pushed apart by steric repulsion, the angles formed by hydrogen bonded atoms becomes smaller [ Fig. 1(C,D)]. This effect originates from the protein structure-derived hydrogen bond potential for helices in Rosetta. Sidechain hydrogen bond geometries in protein structures agree closely with those expected from MP2 QM energy landscape calculations, but helix hydrogen bonds differ considerably because of the constraints presented by protein backbone geometry. 29 The error is that the protein structurederived potential already incorporates the effects of steric and other interactions, so that including both the helical structure-derived hydrogen bond potential and the Lennard-Jones potential results in double counting.
We correct for this double counting by subtracting the contribution of the Lennard-Jones interactions from the angular dependence of the helix backbone hydrogen bond using the iterative procedure described in the Methods. The corrected potential is now more favorable for NÀ ÀHÀ ÀO angles near 1808 and HÀ ÀOÀ ÀC angle near 1658. An additional ensemble of structures was generated with the corrected potential, and the NÀ ÀHÀ ÀO angle peak was found to be 1628 and the HÀ ÀOÀ ÀC angle peak at 1508, matching crystal structures. In addition, the peak OÀ ÀO distance shifted to 3.4 Å , agreeing with crystal structures as well. This improved agreement in a feature not directly controlled by the modification of the potential suggests that the correction strategy is on the right track. There are remaining errors in the region HÀ ÀOÀ ÀC > 1608 [ Fig. 1(D)] with the corrected potential, because corrections are applied only to the position of the peak; the remaining difference in the distribution is in regions with very low populations.

b-sheet CaÀ ÀO hydrogen bond
The first peak in the backbone b-sheet NÀ ÀO or NÀ ÀN distance distribution is very similar in crystal structures and low-energy rosetta models. In addition, the backbone b-sheet hydrogen bond angle distributions in Rosetta models match quite closely to those in crystal structures. Thus, the geometry immediately around the  b-sheet hydrogen bond appears to be well modeled in Rosetta. However, there is a discrepancy in the atom pair distributions in b-sheet structures that is evident in the second peak of the NÀ ÀO and CaÀ ÀO distributions. The CaÀ ÀO distribution shows a peak at 3.3 Å in the crystal structures (Fig. 2), which is somewhat flattened and pushed further away in Rosetta models. The NÀ ÀO distribution in crystal structures has a clear peak between 4 and 5 Å [ Fig. 2(A)], whereas Rosetta models have a more flattened distribution.
We considered the possible origins of this discrepancy. There have been many suggestions that nonpolar hydrogen bond interactions contribute to protein structural features. 30-32 In antiparallel b-strands, the shearing of the neighboring strands leads to nonlinear COÀ ÀNH hydrogen bonds. 33,34 This was proposed to be due to  the nonpolar hydrogen bond interactions. 35 The propensity of right hand twist in b-sheets has also been suggested to arise from Ca hydrogen bond interactions. 30,35 The peak in the CaÀ ÀO distribution in crystal structures could reflect contributions of such hydrogen bonds, or it could be a secondary consequence of b-sheet packing and other well-understood contributions to protein energetics. However, the Rosetta forcefield explicitly models b-strand hydrogen bonding and accounts for the major contributions to protein energetics, hence the former explanation seemed more likely.
We tested the incorporation into the forcefield of a simple distance-dependent nonpolar hydrogen bond potential between the backbone carbonyl oxygen and Ha. The starting guess at the form of the potential was based on the logarithm of the distance distribution between the pair of atoms in crystal structures, normalized such that the potential is 20.5 at 3.0 Å and goes to zero at 3.6 Å . This starting guess was then refined by application of the iterative correction protocol described in the methods. The converged OÀ ÀHa potential has a depth about 40% the strength of the polar hydrogen bond.
With the incorporation of the nonpolar hydrogen bond term, the b-sheet geometry of low-energy Rosetta models closely matches crystal structures. The radial distribution of O to Ca now has the clear first peak found in crystal structures, and more significantly (since it is not explicitly enforced by the new term), the peak in the NÀ ÀO distribution between 4 and 5 Å is now also clearly distinguishable. The registration of paired b-sheets also becomes more similar to crystal structures. As shown in Figure 2(E), the distribution of r, the angle between strand dimers in adjacent sheets 36 used in the Rosetta low-resolution energy function, is shifted relative to crystal structures with the standard Rosetta energy function. With the incorporation of the nonpolar hydrogen bond, b-sheet registration in Rosetta models closely matches that in crystal structures. This agreement is again significant, because the r distribution improvement is not forced directly by the nonpolar hydrogen bond term; instead the improvement in registration and the NÀ ÀO distribution suggests b-sheet geometry as a whole has improved, and further that the nonpolar hydrogen bond contributes to the observed geometry in protein structures.
This example illustrates a strength of our approachwe make inferences not from the distributions observed in crystal structures, as for example traditional knowledge-based potentials-but from the differences between these distributions and those of models generated with a forcefield representing the major contributions to protein energetics. Differences as in this case are a strong argument for missing physical chemistry.

Backbone torsion potential
For each amino acid in each secondary structure (helix, strand, turn), we collected statistics on the Rama-chandran (u,w) backbone torsion angle distribution in crystal structures and low-energy Rosetta models and determined those for which the two were most in disagreement. The largest differences were found in the w angle distribution in b-sheets (Supporting Information Fig. S1). In crystal structures, the w angle distribution in b-sheets for most residues has a single peak centered around 1408; lower values of w are disfavored due to steric repulsion between Cb and O (Fig. 3). However, for polar and charged residues with short sidechains, Asp, Asn, Thr, and Cys, there are two distinguishable peaks in the w distribution in b-sheets, one around 1208 and the other around 1408$ 1508. The peak around 1208 arises because the steric repulsion between Cb and O is compensated by a hydrogen bond formed between the sidechain and the backbone amide of residue i 1 2, which constrains the backbone w angle to be lower than 1208 [ Fig. 3(B)]. The correlation between this hydrogen bond and w is evident in crystal structures, especially in the loop region as shown in Figure 4. A peak in distribution arises at w around 1208 and distance below 3 Å between Asn Od 1 and N of residue i 1 2.
Rosetta models the w distribution well for most residues (Supporting Information Fig. S1), reproducing the single peak around 1408 and matching the intensity. However for Asn and Asp, Rosetta underestimates the population around 1408, while oversampling the region with w < 1208. The discrepancy between Rosetta models and crystal structures is due to energy overcounting. In the Rosetta forcefield, this conformation is favored by three energy terms, the hydrogen bond potential, the Ramachandran potential, and the backbone-dependent rotamer potential. The favorable Ramachandran and rotamer potentials both reflect the frequent occurrence of the hydrogen bond in crystal structures, and this double counting results in overstabilization of the w < 1208 conformation. A clear artifact is that the Asn conformation is still favored even when the hydrogen bond is not formed. This leads to a correlation between the atom pair distance between Asn Od1 and N of residue i 1 2 at over 4 Å and the w angle < 1208 in b-sheets [ Fig. 4(D)]. That the origin of this artifact is the overcounting is confirmed by the observation that when Ramachandran potentials of Asn and Asp are substituted with those of Ala, the peak of w below 1208 disappears.
We corrected for this double counting by iteratively correcting both the Ramachandran potential and the rotamer potential as described in the Methods. As two terms are being modified, there could be multiple independent ways to achieve the same net correction-to resolve this degeneracy, we aimed for corrections which made both the backbone torsion potential and the rotamer potential more similar to those of other residues. With the corrections, the backbone torsion potential for Asn becomes more similar to those of nonpolar residues, with the energy in the w < 1208 region less favorable [ Fig. 3(E)]. The rotamer potential also no longer favors the v 2 $ 158 region when w is less than 1208, becoming less backbone dependent. With the optimized potential, the sidechain-backbone hydrogen bonds in the w < 1208 still form, but they are only favored by the hydrogen bonding potential. The backbone torsion angle distributions in the modeled structures using the optimized potentials are more similar to the crystal structure distributions (Figs. 3 and 4) than those using the original potentials.

Coupling between Lennard-Jones interactions and sidechain torsion potential
Comparison of sidechain rotamer distributions between low-energy Rosetta models and crystal structures shows that the distributions of rotamers at v 1 $ 2608 are underpopulated in Rosetta models for u < 21208  2608). Structural analysis shows that as u shifts below 21408, the backbone carbonyl group of residue (i 2 1) moves toward the sidechain, leading to steric repulsion. This clash is small, so the repulsion does not completely disfavor rotamers with v 1 $ 2608, as hey  show that Rosetta stabilizes the conformation in Fig. 3B as in crystal structures, where the hydrogen bond restrains the backbone w to lower than 1208. The blue arrow in (D) highlights the artefact with the standard Rosetta forcefield; the hydrogen bond is not formed, yet the backbone w is still favored to be lower than 1208. The sampling of this region is much reduced with the optimized forcefield in (F). The color scale shows the minus log of probability of a given w and Od1À ÀN distance bin, offset by E mean [Eq. (6)].
are still observed in crystal structures, but this sidechain rotamer is significantly more depleted in low-energy Rosetta models (the rotamer is observed in 20$30% of crystal structure positions with a u < 21408, whereas in low-energy Rosetta structures, the frequency is <5%).
Double counting energy contributions are again likely to be the origin of the discrepancy: both the rotamer potential and the Lennard-Jones potential disfavor rotamers with v 1 $ 2608, when u is less than 21408. To remove the double counting, the rotamer potential was optimized using the iterative protocol described in the Methods section to effectively subtract the Lennard-Jones contribution from the rotamer potential. As shown in Figure 5, the energy at u < 21408 is unfavorable in the original rotamer potential. After optimization, for the same w, the potential at u < 21408 and u > 21408 is now similar for the v 1 $ 2608 rotamer [ Fig. 5(B)]. With the corrected rotamer potential, the v 1 $ 2608 rotamers of Val are now properly populated; the distribution in low-energy Rosetta models is now much closer to the crystal structure distribution (Fig. 5).
Additional errors were found in Tyr and Phe for which v 2 near 0 is overpopulated in the low-energy Rosetta mod-els. For example, for buried positions, the Phe rotamer at v 1 5 -708 and v 2 5 -148 has an occupancy of 8% in crystal structures but over 20% in Rosetta models. The origin of the discrepancy here is more difficult to identify because the rotamer distribution varies little with all Rosetta energy terms up-or down-weighted individually. However, for exposed Phe with more than 20% solvent exposure, this rotamer is 9% occupied in Rosetta models, matching the rotamer library well. This suggests that packing of the rest of the protein suppresses rotamers with v 2 near 908 and enhances the probability of rotamers with v 2 near 08. After optimizing the rotamer potential using the iterative protocol, the energy of the v 2 near 0 rotamers increases significantly so that the correct rotamer distribution is now recovered in Rosetta models.

Impact of potential function corrections on the overall energy landscape
One of the key challenges in forcefield optimization is to prevent over-optimizing toward one set of measurements, while worsening those on others. We tested the  combined changes described above by generating with the new forcefield extensive energy landscapes for each of the 110 proteins in our test set. The energy gap between near-native structures and low-energy computed structures was calculated as described in the Methods section. For each protein, the energy gap with the optimized forcefield is compared to that with the standard forcefield in Figure 6. The overall impact on the energy landscape is small, on average the energy gap is shifted by 0.5 Rosetta energy unit (1 Rosetta energy unit $ 0.5 Kcal/mol) in favor of native structures. Of 110 proteins tested here, in 19 cases, the changes in the energy gaps now favor the native structure more by 1.5 Rosetta units. No protein has an energy gap changed by more than 3 energy units. Thus, the optimization of the forcefield does not significantly alter the energy difference between native and non-native structures, while improving the geometry.

Independent benchmark
The optimized Rosetta forcefield was applied to an additional independent benchmark, using 55 protein structures from the CASP8 experiment 25 as described in the Methods section [ Fig. 6(B)]. The optimized forcefield has a small impact on the energy gaps of these proteins overall. On average the native structures are more favored by 1.4 Rosetta energy units. There are five proteins that have energy gaps > 210 Rosetta units, where the discrimination between native and decoy structure is poor. All of these proteins show improvements over 2 Rosetta units favoring the native structure.

DISCUSSION
We show how comparison between distributions of structural features in crystal structures and low-energy computed structures can be used to guide forcefield optimization. Errors in the forcefield are detected in regions where the distributions in low-energy models differ from those in crystal structures. We apply the method to resolving the interdependencies between sidechain and backbone torsion potentials and backbone and sidechain hydrogen bonding interactions, and the approach also motivates explicit treatment of Ca hydrogen bonds. The new forcefield yields minima with improved geometry without significantly changing the position or relative depth of minima on the overall energy landscape. The approach is physically based and structure guided and can be applied more generally to improve forcefields using information from macromolecular structures.
It is instructive to compare the approach to correcting torsion potentials described here to the derivation of torsion potentials from quantum chemistry energy calculations on small peptides in MM forcefields. In both Rosetta and MM forcefields, the torsion potential is something of a catch-all for energy contributions not contained within the other more easily parameterizable terms in the potential. MM methods parameterize the torsion potential by computing the total energy of small peptides for varying values of the torsion angles and then subtract the contribution of all other terms-the remainder is the torsion potential. Likewise, we optimize the torsion potential to reproduce the correct distribution keeping other terms fixed. Both approaches have strengths and weaknesses; the QM approach has the advantage of parameterizing on energies rather than frequencies, but the disadvantage that the energy calcula- Comparison of energy gap between the native and non-native structures with optimized forcefield versus standard Rosetta forcefield (1 Rosetta energy unit % 0.5 Kcal/mol). A more negative energy gap corresponds to a better discrimation between the native and decoy structures. (A) 110 proteins from the dataset that the forcefiled optimization is carried out on and (B) 55 proteins in the CASP8 benchmark set. For each protein, 10,000 full atom structures were generated with Rosetta refinement protocol using the standard or optimized Rosetta forcefield, and energy gaps between the lowest native and non-native conformations were compared as described in the Methods section. For each protein in the sets, the energy gap with the optimized forcefield is plotted against the energy gap with the standard forcefield. The dashd lines indicate that energy gap change by 1.5 energy units (A) or 5 energy units (B). Although the energy gaps are quite similar, for most difficult discrimination problems (energy gaps close to 0), there is small but consistent improvent with the optimized forcefield. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] tions can only be done on very short peptides and not in explicit solvent, while the distributions are collected from proteins of all sizes in real water. A further parallel is the adjustment of the backbone torsion potential in going from AMBER94 to AMBER99SB to favor b-strand forming peptides 13 and sidechain torsion potential improvement in AMBER99SB-ILDN 37 ; our approach has the advantage of using samples spread throughout the energy landscapes of large numbers of proteins.