Improved prediction of protein side-chain conformations with SCWRL4

Input and construction of coordinates

An individual residue position is defined by specifying four backbone atoms (N, Cα, C, O) in a PDB-format input file. These individual residue sites can comprise one or more polypeptide chains, from which the backbone dihedral angles ϕ and ψ are calculated for each residue. For purposes of looking up residues in the rotamer library, the N-terminal residue ϕ is set to −60°. Similarly for C-terminal residues, ψ is set to 60°. These values are those for which there is weak dependence of the rotamer probabilities on the missing dihedral angle.38 The C_i−1N_i atom distances are checked to determine whether there are missing internal residues in a chain.

For each residue, rotamers are read from a new version of the backbone-dependent rotamer library (Shapovalov and Dunbrack, in preparation). This rotamer library is based on a much larger data set, and is derived using kernel density estimates and kernel regressions. The rotamer library includes rotamer frequencies and mean dihedral angles and their standard deviations over a discrete (ϕ, ψ)-grid. This library offers much greater detail for nonrotameric degrees of freedom, in particular χ₂ for Asn, Asp, His, Phe, Trp, and Tyr and χ₃ for Glu and Gln. Optionally SCWRL4 can determine frequencies and dihedral angle parameters by bilinear interpolation from the four neighboring ϕ, ψ grid points in the library. For each χ₁ rotamer of Ser and Thr, SCWRL4 generates three rotamers for the hydroxyl hydrogen with χ₂ dihedral set to −60°, +60°, and 180° and the variance set to 10° times the corresponding parameter given in the configuration file. For each χ_1, χ₂ rotamer of Tyr, two rotamers are generated for the hydroxyl hydrogen with χ₃ dihedral set to 0 and 180°, which are the values observed in neutron diffraction studies.39 For His, extra rotamers are created for the singly protonated states (proton on ND1 or NE2). Rotamers that represent positively charged His can be enabled in the program using an option in the configuration file.

Side-chain coordinates are built for all rotamers and for subrotamers about these rotamers used by the flexible rotamer model (FRM, see later). Subrotamers as used here are conformations with dihedral angles ± one standard deviation (or a fixed proportion thereof) away from rotamer values given in the rotamer library. For subrotamers, only one dihedral at a time differs from the library value, since we found that allowing multiple deviations did not noticeably improve the accuracy but did slow the calculation (data not shown). Side chains are represented in a tree-like structure, so that atoms common to more than one subrotamer (e.g., same CG position for different χ₂ conformers) are calculated and stored only once.40 Coordinates are built using a fast incremental torsion to Cartesian conversion method.41

Because SCWRL4 uses a large number of rotamers and subrotamers, we implemented a fast collision detection algorithm based on k-dimensional Discrete Oriented Polytopes or kDOPs.37 The kDOP algorithm is based on two key ideas. The first idea is to enclose each geometric object into a convex polytope of a special kind and use these as bounding boxes for clash checks. A particular class of kDOPs is defined by a set of k pairwise noncollinear unit vectors, and consists of all convex polytopes with 2k facets such that any facet is perpendicular to one of these vectors. Examples are shown in Figure 2. For instance, if k = 3, these could be the x, y, and z-axes, so that all bounding boxes are rectangular parallelepipeds whose faces are perpendicular to the Cartesian axes. The second key idea is to organize all bounding boxes related to a particular group of geometric objects as a hierarchy. These hierarchies can then be used for efficient search of all possible clashes between individual bounding boxes.

The advantage of using a single set of vectors for defining the bounding boxes is that two bounding boxes of the same kDOP class can be efficiently checked for clashes. As illustrated in Figure 2, this can be accomplished by testing for overlaps of the real intervals that represent their projections across the corresponding basic axes. Let {a_i} and {b_i} be the sets of these intervals for two kDOPs “A” and “B” respectively. “A” does not clash with “B” if there exists i ∈ {1..k} such that min (a_i) > max (b_i) or if there exists j ∈ {1..k} such that min (b_j) > max (a_j). If neither of these conditions are met, then the underlying objects enclosed inside of these two boxes may clash, but do not necessarily do so. In the case of atoms, this is to be checked by pairwise distance calculations of the objects enclosed in the bounding boxes.

Building a kDOP around a geometric object consists of finding projections of the object onto the basic axes. Both the van der Waals function and the hydrogen bond potentials described in the next section have a certain boundary distance beyond which the potential is zero. These distances are used to represent each atom as a sphere with a certain radius. To build a bounding box around a whole side chain, each atom is enclosed into a kDOP and then the elementary shells are merged. The rotamers and subrotamers of a side chain can be enclosed into a single kDOP, such that all residue–residue interactions can be checked very quickly. In SCWRL4, we use four basic axes and construct bounding boxes around individual atoms, backbone atoms of each residue, parts of each side chain, individual rotamers (i.e. entire set of its subrotamers) and every residue (all rotamers). The basic axes form a tetrahedral geometry:

Using these four axes results in somewhat faster calculations, by about 15%, than using three axes along the x-, y-, and z-axes, despite some overhead involved in calculating the projections.

Calculation of self-energies and pairwise energies via modified FRM

SCWRL4 uses both a rigid rotamer model (RRM), as in SCWRL3, and a FRM.32 In the RRM, the total energy of the system is expressed as:

where the vector r specifies a single rotamer for each of N residues in the system. In this case, the self-energy of each rotamer is:

where the first term expresses the rotamer energy relative to the most populated rotamer, r_max, given the backbone dihedrals ϕ and ψ of residue i and the frame term expresses interaction of the side chain with the backbone and any ligand or other fixed atoms present. We allow the value of the constant in front of the log term to be residue-type dependent.

In contrast to SCWRL3, in SCWRL4 the frame and pairwise rotamer energies consist of repulsive and attractive van der Waals terms as well as a hydrogen bonding term. The repulsive van der Waals term is the same as the piecewise linear term used in SCWRL3, but is combined with a short-range attractive potential as follows. If σ_ij is the sum of the hard-sphere radii of atoms i and j and E_ij is , where the E_i values are the E_min values from the CHARMM param19 potential,42 and d is the distance between the two atoms, then

This potential is shown in Figure 3 along with the standard Lennard–Jones potential with the same E_ij and R_ij. The hard-sphere radii were manually optimized for the training set accuracies. The minimum energy occurs at which is close to the standard Lennard–Jones parameterization in which the minimum occurs at r_min = 2^1/6σ_ij = 1.12σ_ij. The parameters are given in Supporting Information.

The hydrogen bonding term in SCWRL4 is similar to the one used in Rosetta,36 although it is parameterized in a different way, as shown in Figure 4. We define d in this case as the distance between a polar hydrogen (HN- or HO-) and a hydrogen bond acceptor (oxygen), as a unit vector from O acceptor to H, as a unit vector along the covalent bond from the hydrogen bond donor heavy atom D to H, and two unit vectors and from the hydrogen bond acceptor O toward the middle of the oxygen lone-pair electron clouds. For carbonyl oxygen, these two vectors are 120° apart from the double bond and coplanar with the carbonyl carbon substituents. For hydroxyl oxygen, these two vectors are 109.5° from each other and from the other two oxygen substituents (H and C), forming a tetrahedral arrangement. The hydrogen bond function is evaluated first for , and if no hydrogen bond is found, then for . For , the weight w for the hydrogen bond is defined as:

where is the angle between the DH bond and the H…O vector and is the angle between the O-lone pair and the O…H vector. If the multiplicand under the square root is negative then the score is set to zero. The calculation of this score enables an efficient implementation and together with the distance d and vector can be done within 30 arithmetic operations, one division, and two square root evaluations.

After the weight w has been computed, it is used to derive the final energy of oxygen-hydrogen interaction by balancing the default van der Waals energy and pure hydrogen-bond attraction terms:

where q_H and q_O are the charges from the CHARMM param19 potential. The formulas mentioned earlier include five atom-independent coefficients: d₀, σ_d, α_max, β_max, B. The values of these coefficients were optimized on the training set proteins and are given in Supporting Information.

Using single rotamers sometimes results in poor packing predictions, due to fluctuations in the dihedral angles and imprecise representations of the backbone in homology modeling. We investigated the use of subrotamers, which we define as conformations that differ in one or more dihedral angles by one standard deviation (or some constant times this value) from the mean values given in the rotamer library:

If we allowed variations in all dihedral angles in this manner, treating the subrotamers as additional rotamers resulted in intractable calculations using the graph decomposition algorithm used in SCWRL3. Even with the tree decomposition of Xu,31 implemented in SCWRL4 (see later), the calculations often remained intractable. So we implemented the FRM of Mendes et al.,32 in which the subrotamers are integrated to produce an approximate free energy using the Kirkwood superposition approximation43:

We treat the first term as the “self-free energy” and the second term as the “pairwise free energy.” A_self(r_i) and A_pair(r_i,r_j) are defined as:

The terms E_frame(r_i,s_i) and A_frame(r_i) contain only the van der Waals and hydrogen bond energies. In our implementation, each residue type has a separately optimized temperature, and for the pairwise free energy, T_ij = (T_i + T_j)/2.

Graph construction

As with SCWRL3, some rotamers with high self-energy are removed from the calculation, since they are very unlikely to be part of the predicted structure. These rotamers are marked as inactive. In SCWRL3, rotamers with self-energy above a certain residue-independent bound were inactivated. However, it sometimes happens that all rotamers have self-energy above this bound. In this case all rotamers were reactivated. In SCWRL4, we replaced this heuristic by making the bound relative to the lowest energy rotamer for each residue. This approach guarantees that at least one rotamer will remain active. After some study the value of this threshold was set to 30. The exact value of the threshold can be customized through the configuration file.

Before the graph is constructed, disulfide bonds are resolved. SCWRL4 uses the same criterion as SCWRL3 to identify if two cysteine side chains can form a disulfide bond, but introduces a new procedure for resolving ambiguities. An ambiguity occurs when more than one rotamer of a particular cysteine residue can form a disulfide bond or when one rotamer can form disulfide bonds with more than one other cysteine side chain. To select a particular collision-free combination of disulfide bonds, SCWRL4 finds the minimum total energy out of all possible combinations of feasible disulfide bonds. To do this, we use an objective function in the form:

where the summations run over all possible disulfide bonds, C is a large positive constant and η is a binary function that evaluates to one if a particular disulfide bond is switched on and to zero otherwise. The functional above is of the same form as the one used to compute the total energy of rotamer assignment. Therefore we can minimize it for function η via the same optimization procedure. Doing this yields, a list of the optimal disulfide bonds that do not have collisions. If for a particular cysteine residue one of the rotamers is part of an optimal disulfide bond then all other rotamers are inactivated for that residue. Energies of interaction of cysteines in disulfides are added to the self-energies of rotamers of other side chains within interacting distance.

SCWRL uses an interaction graph to represent the side-chain placement problem.17, 30 In this graph, vertices represent residues while edges between vertices indicate that at least one rotamer of one residue has a nonzero interaction with rotamers from another residue connected by the edge. For a single protein or protein complex, the graph is constructed by checking for overlap of the kDOP around whole residues. If at least one rotamer or subrotamer of one residue can interact with nonzero energy with a rotamer or subrotamer of another residue, then an edge is added to the graph between the vertices in the graph representing these residues.

SCWRL4 is able to model side chains in symmetric complexes using symmetry operators. These rotation-translation operators can be generated from the CRYST1 record in the input PDB file or specified explicitly by the user in a separate input file. For crystals, if the input PDB file contains the asymmetric unit, all residues in asymmetric units that may contact the input coordinates are constructed, as described previously.44 Bounding boxes are constructed around the residues, rotamers, subrotamers, and atoms of the symmetry copies. Interactions between atoms in the input structure and its side chains and atoms in the symmetry copies and their side chains are determined. If side chains in the input structure interact with the backbone or ligands of the symmetry copies, then the static frame energies of these residues are modified accordingly. If the side chains of the input structure interact with the side chains of the symmetry copies, then an edge is created between the corresponding residues, if it does not already exist. If it does, then the pairwise energies are modified to account for the additional interactions between symmetry-related proteins. Thus a residue on one side of a protein may have an edge with a residue on the other side of a protein because of symmetry.

Graph solution via tree decomposition

Before the major optimization via dynamic programming is launched the interaction graph undergoes some preprocessing consisting of edge decomposition and DEE. Typically, this eliminates a significant number of rotamers as well as some edges and nodes. Because edges were formed based on overlapping of bounding boxes some of them may contain only zeros as pairwise rotamer–rotamer energies. If this is the case or if the actual energies of interactions are very close to zero then the edge is removed.

Edge decomposition removes edges that can be approximated as the sum over single-residue energies. If this representation is feasible within a certain threshold then the corresponding self-energies are modified and the edge is removed. With larger thresholds, more edges may be removed. In this preprocessing stage, the threshold is set to a very small value, ϵ = 0.02 kcal/mol.

The pairwise energies of two residues, E_pair(r_i,r_j), in the RRM or free energies, A_pair(r_i,r_j), in the FRM, may be represented by a matrix of real numbers e_kl for rotamers k = 1…m and l = 1…n. Edge decomposition consists of finding two sets of real numbers {a_k} and {b_l} which minimize the average deviation:

By setting the partial derivatives of δ with respect to a_k and b_l to zero, we find that these two sets should satisfy the following equations:

The initial task is not well defined as its solution is not unique. Thus adding some value to all a_k and subtracting the same value from all b_l does not change the sum a_k + b_l. Therefore, we can set to an arbitrary value. For example, we can set:

Substituting this value into the second equation, we find the corresponding value for :

Using these values, we can determine a_k and b_l and evaluate the maximal absolute deviation:

SCWRL4 checks if this deviation is less than a certain threshold and if so then it removes the corresponding edge and modifies the self-energies of the kth rotamer of residue i and the lth rotamer of residue j:

As stated earlier, the initial value of the threshold is 0.02, which enables the algorithm to eliminate almost all redundant “near-zero” edges. We remove from the graph those nodes that now have zero edges; its assigned rotamer is that of lowest E_self.

The next step is to perform DEE that identifies and removes rotamers that cannot be the part of the global solution. These rotamers are identified via Goldstein's criterion that was used in SCWRL3.27 If for a certain residue only one rotamer is left after DEE then that rotamer is part of the solution. If this residue has adjacent edges then all pairwise energies with the remaining rotamer are incorporated into self-energies of the corresponding rotamers from adjacent residues and these edges are removed. This makes the residue isolated, which means that it can be removed from the graph; the self-energy of its single rotamer is added to the total value of the minimal energy. The edge decomposition and DEE steps are repeated until nothing further is removed.

As in SCWRL3, the resulting graph may contain separated subgraphs or clusters with no edges between them; each of these clusters is then subject to graph decomposition. In SCWRL3, the graph decomposition was based on the determination of biconnected components, which are subgraphs that cannot be broken into parts by the removal of a single node. The graph is then a set of biconnected components connected by single nodes called articulation points. Tree decomposition can be viewed as a generalization of graph decomposition based on biconnected components.31 To see this, we show in Figure 5 the same graph as described in the SCWRL3 paper, its decomposition into biconnected components, and its tree decomposition. The nodes of the graph on the left are gathered into “bags” which are nodes of the tree shown on the right. Every node of the graph is represented in one or more of the bags (definition of Condition 1 given later). Every edge of the graph on the left is also represented in one or more of the bags, so that the two nodes of an edge are together in at least one bag (Condition 2). Finally, for any vertex of the graph, all those bags on the tree that contain the vertex form a connected subtree (Condition 3). More formally:

Definition

A tree-decomposition of a graph G = (V,E) is a pair (T,Z), where T = (W,F) is a tree (i.e., a graph with no cycles) with vertex set W and edge set F and Z = {Z_w : Z_w ⊆ V}_w∈W is a family of subsets of the set V associated one-to-one with the vertices of T that satisfies the conditions:

Because of the one-to-one correspondence between sets Z and W, we will denote the vertices of tree T as “bags.” Figure 5 shows that Condition 1 is satisfied by this tree decomposition, since all the residues are present in one or more bags. Residues c,d illustrate that Condition 2 is satisfied since the edge c-d is contained in at least one bag (in this case, two). Residue h illustrates Condition 3, since all the bags that contain h are connected in a single subtree.

Typically several different tree-decompositions can be built for a given graph. The width of a particular tree-decomposition is the size of the largest bag minus one. For a given graph a tree-decomposition with the minimal possible width is the optimal one and its width is called the treewidth of the graph. This characteristic indicates how well a graph is tree-decomposable. For example if a graph has no cycles (and thus is a tree) then its treewidth equals one, while for a simple cycle the treewidth equals two.

In SCWRL3, the graph solution begins with any biconnected component with a single articulation point by finding the minimum energy of the biconnected component residues for each rotamer of the articulation point. This energy is then added to the self-energy of the articulation point rotamer, and the rotamers of the biconnected component that achieve this minimum energy are assigned to the articulation point rotamer. The biconnected component can then be removed, and the process continues for all biconnected components with one articulation point in the remaining graph. The combinatorial problem is thus reduced to the order of the largest biconnected component (i.e., the one with the largest number of rotamer combinations).

In a tree decomposition, instead of using single nodes to separate the graph, the graph can be separated by removing one, two, or more nodes. To see this, in the tree decomposition in Figure 5, each bag w is broken up into two sets of residues, L_w and R_w, where the residues in L_w are those residues in the bag that are shared between the bag and its immediate parent bag. For each bag in the figure, these are listed to the left of a vertical bar. The remaining residues in the bag, the set R_w, are those not in the parent and are placed to the right of the vertical bar. Each set L_w is a separation set of the graph G45; that is removing the residues in L_w breaks the original residue graph into two or more separate unconnected graphs. For instance, removing residues b and c breaks the original graph into two graphs, one consisting of residue a and the other the rest of the graph below residues b and c.

Solving for the minimum energy of the graph proceeds as it does in SCWRL3. Starting at a leaf (a bag with no children), y, e.g. the one containing “b c | a”, find the lowest energy of the residue(s) in R_y (in this case residue a) for each combination of the rotamers in L_y (in this case, residues b and c), saving the corresponding assignment of rotamers of R_y. Then, add these energies to that rotamer combination in the parent bag, which by the definition of tree decomposition contains b and c. The procedure continues up the tree to the parent node of y (let us call it node z). Again, the minimum energy of all the rotamer combinations of those residues in R_z is calculated for each rotamer combination in set L_z. These energies need to include the energies for b,c calculated for the child node y. This procedure continues until only the root bag is left. We provide a more formal description of this procedure later.

The complexity of the solution is associated with the width of the tree, since all the rotamer combinations of the residues in each bag need to be enumerated. It is in general difficult to compute a treewidth and to find the optimal tree-decomposition, and it has been proved to be NP-hard for an arbitrary graph.46 For building a tree-decomposition, we have developed a heuristic algorithm. Our algorithm is similar to the one suggested by Xu31 who referred to it as a “minimal degree heuristic.”

In the first step, the family of sets Z is built. The input graph is gradually disassembled using a loop of the following steps:

• 1

  Select any vertex with the minimal number of adjacent edges.

• 2

  Form a bag of the tree from this vertex and all its neighbors. The selected vertex we will denote as the primary vertex of the corresponding bag.

• 3

  Add edges into the graph being processed so that the neighbors of the selected vertex become a clique (a subgraph where all nodes have edges to each other).

• 4

  Remove the selected vertex and all adjacent edges from the graph.

• 5

  Repeat from the first step until there are no more vertices left in the graph.

Thus, we obtain a set of “bags” which represent the vertices of the tree-decomposition.

Bags are numbered in the order of their construction. It is important to notice here that within any iteration the intersection of the bag w with the vertices of the remaining graph (S_w) consists solely of the neighbors (N_w) of the initial vertex of the bag concerned, Z_w ∩ S_w = N_w. The one-to-one correspondence between vertices and the bags verifies that the first condition in the definition of tree decomposition is automatically satisfied. Also, it is clear that the edges are removed solely during the bag construction and that when any edge is removed both adjacent vertices are included into a bag. This guarantees that the second condition in the definition of a tree-decomposition is satisfied.

The second step is to connect the “bags” to obtain a tree that meets the definition of tree decomposition. This is done by sequentially fastening these bags to the tree in the reverse order in which they were constructed. Thus, the bag that was created last becomes the root of the tree. The next bag becomes connected to it and thus becomes its immediate child. For the next bag, there are two choices for where to attach it. However, the appropriate node of the tree-decomposition must meet the following condition:

According to this condition, a set of vertices in the appropriate node must contain all vertices from the bag to be added that are already present in the tree.

The tree decomposition just obtained undergoes some additional minor processing. This consists of two normalization rules, which are applied until they cannot be applied further: (1) If all vertices associated with some bag belong to the vertex set of its immediate parent then this node is removed and all its immediate children are reconnected to the parent node. (2) If some bag contains all vertices associated with its parent node then the parent bag is substituted by this bag which thus moves up one edge towards the root.

The minimum energy rotamer configuration is calculated as follows. Starting with a leaf node consisting of sets L and R, the left and right portions of each bag as defined earlier, we define 〈L〉 as the set of all rotamer combinations of the residues in the set L, and similarly define 〈R〉 for set R. A single member of 〈L〉 we denote as l, which is a vector of rotamer assignments, one rotamer l_i for each residue i in the set L; similarly define r for 〈R〉. For a leaf node, we calculate energies for each vector l:

where

and

For fixed rotamers in L, only the pairwise interactions with rotamers in R are included, while both self and pairwise interactions among the rotamers in R must be added. For the FRM, the values of A_self and A_pair are used instead of E_self and E_pair. For an inner node of the tree decomposition, we need to add in the energies assigned to rotamer combinations in the node via its children:

where

where the sum is over the children c of the inner node, and l_c is a vector of the rotamers of the residues in the set L_c, such that these rotamer are in the set {l_i : i ∈ L; r_j : j ∈ R}. By definition of the tree decomposition, all the residues in L_c are in {L ∪ R}. To calculate the energies of the internal nodes, the nodes must be traversed in leaf to root order. Since the root has no parent, it has no set L (equivalently, its set L is empty). Its energies are given by

In the last part of the algorithm, the nodes of the tree-decomposition are traversed in the root-to-leaves order to assemble the global assignment of rotamers for the cluster being processed. For any node except the root, we have a local partial solution that lets us obtain an optimal assignment of rotamers of R for any rotamer assignment of rotamers over L. But by construction L belongs to the parent bag, which means that we can easily retrieve the optimal assignment over all residues in some bag if we know the optimal assignment at the parent node. Thus, we have a recursive procedure that gradually extends an assignment for the entire cluster starting from the root node:

For each child c of the root, the optimal assignment of rotamers to the residues in R_c may be made, given the assignment of rotamers of the root, and the minimum energy added to the total:

where the assignments in l_c are already known from r_root. The rotamers are assigned and the energies updated for each child of each c, and so on, following from the root to all leaves in a depth-first search order.

The actual search for the local solution is done via exhaustive direct enumeration, which is quite affordable if the product of rotamer numbers within the corresponding bag is not very large. The number of possible rotamer assignments over a particular bag we will refer to as local complexity of the node. The sum of the local complexities of all nodes gives the overall computational complexity of the optimization. Typically, tree-decompositions of smaller width yield lower complexity. To limit the time required by SCWRL4 for a single rotamer assignment we introduced an upper bound for the overall complexity of 10⁸. If the actual complexity exceeds this bound then the optimization is treated as not tractable and SCWRL4 returns to the graph construction step (edge decomposition and DEE) after doubling the value of the edge decomposition threshold. This process continues until a solution is found.

Output of the results

The optimization resolves both the minimal total energy of the entire model and the corresponding assignment of rotamers. The SCWRL4 executable saves the resolved optimal conformation of the whole protein model into PDB file. The corresponding value of the total energy is printed into the standard output, which can be redirected to a file for further analysis. If the task was set up and solved via the API of the SCWRL4 library then the corresponding workspace with or without modifications can be used in subsequent calculations.

Training and test sets

We constructed a training set of proteins for optimizing the parameters and procedures, and a separate testing set for reporting the accuracy of SCWRL4. Because we wanted to use electron density calculations to estimate the reliability of side-chain coordinates, we started with the list of PDB entries with electron densities available from the Uppsala Electron Density Server,47 generally those with deposited structure factors. We removed entries with ligands other than water, so that side chains could be predicted without requiring charges, hydrogen positions, or van der Waals radii of ligands. This set was culled using the PISCES server48, 49 at maximum mutual sequence identity 30%, ≤1.8 Å resolution, and maximum R-factor of 20%. Because we planned to optimize the energy function by predicting side chains in the crystal form, we checked whether the CRYST1 records and scale matrices produced viable crystals. We removed some entries that produced extensive clashes of protein atoms when crystal neighbors of the asymmetric unit were constructed (e.g. PDB entry 1RWR). The resulting list of proteins was broken up into the training and testing set, with the training set consisting of monomeric asymmetric units to speed the optimization procedures described later.

For all complete side chains in the resulting protein lists, we calculated the geometric mean of the electron density, as described previously.50 In this prior work, low values of mean electron density were correlated with nonrotameric side chains and conformationally disordered side chains. For each residue type, mean electron densities for the training and testing sets were sorted and turned into percentiles, with 0% for the lowest electron density side chain and 100% for the highest.

For both sets, we used the program SIOCS (Heisen and Sheldrick, unpublished) available with SHELX51 to resolve the ambiguity in the flip state of Asn and His χ₂ and Gln χ₃. SIOCS uses hydrogen bonding and crystal contacts to indicate whether these side chains are correctly placed in crystal structures, or if it is likely that the terminal dihedrals should be flipped by 180°.

In crystal mode, calculations were performed on the asymmetric unit with inclusion of interactions with crystal neighbors. However, accuracy was only assessed on one protein of the asymmetric unit, as provided by PISCES.

Optimization of the energy parameters

The accuracy of SCWRL4 depends on the choice of rotamer library, the nonbonded energy functions and their parameters, the parameters used in the FRM, and several other procedural choices. We first need to decide on an objective function to be optimized. There are a number of possible choices, including RMSD, percent χ₁ correct within some threshold (typically 40° in previous literature), percent χ₁₊₂ correct, and percent of side chains correct (side chains with all χ angles within some threshold). We decided to use the average absolute accuracy. For a side-chain type such as Lys, this is an average of percent χ₁, percent χ₁₊₂, percent χ₁₊₂₊₃, and percent χ_1+2+3+4 correct:

where N₁₂ for instance is the number of lysine side chains with both χ₁ and χ₂ correct within 40°. This value gives added weight to the more reliably determined degrees of freedom closer to the backbone. To obtain accuracy across all side-chain types, we weight PC for each amino acid type by its frequency:

We have a large number of parameters that can be optimized. First, the χ-angle deviations δ_i for the subrotamers can be set individually for each dihedral degree of freedom. SCWRL4 uses constant times the standard deviation provided in the rotamer library, where the constant is specific for each degree of freedom for each side-chain type:

We also optimize the “temperature” used in the FRM procedure separately for each amino acid type. For calculation of pairwise rotamer–rotamer energies, the temperature is taken as the arithmetic average of the temperatures of the corresponding amino-acid types. The last parameter is the coefficient in front of the rotamer library term which balances the influence of the static frame and the rotamer library in the self-energy of a rotamer. Thus for every amino-acid type, we obtain l + 2 parameters, where l is the number of χ angles required to specify the conformation of the side-chain of a certain amino-acid type. These parameters form a 78-dimensional space that can be searched to improve the quality of prediction. In addition, we also have the hydrogen bond parameters and the atomic radii that can be optimized to improve the prediction accuracy.

We used a special technique to perform the optimization in the 78-dimensional space of the parameters concerned. The main idea is similar to classical block-coordinate descent methods,52 where on each iteration an optimum along one or several axes is resolved. In our case, the search space of any iteration consists of the parameters of a single amino-acid type, while the objective function is maximized within some vicinity of the current values of these parameters. Every iteration update only the parameters associated with the corresponding amino-acid type whereas the parameters of the other amino-acid types remain unchanged. In this method, any iteration updates the parameters in a way that increases the value of the objective functions, otherwise keeping the parameters unchanged. Within each iteration, the approximate solution of the underlying optimization task is obtained using a special technique as described in Supporting Information. Changing the amino-acid types from one iteration to another will impose a sequence of points in the original 78-dimensional space along which the value of the objective function will gradually increase (or at least will remain unchanged). Iterations are grouped into rounds—a randomly shuffled sequence of 18 iterations, which contains all residue types except ALA and GLY. Our experiments showed that a few rounds are typically enough to find some (not necessarily global) maximum of the objective function.

The key advantage of the protocol is that it lets us use specific properties of the internal structure of the rotamer assignment procedure. When the FRM is enabled, the most CPU-consuming part of the whole prediction is the calculation of the interaction graph and especially the calculation of the pairwise energies. However for some pair of amino acids, if neither of the rotamers has been changed then the current energy of the pairwise interaction between them is valid and does not need to be recalculated. Conversely, changing the parameters of some amino acid type affects only the rotamers of that type and their interactions with rotamers of any other type. During a single iteration, only side chains of one amino-acid type are modified and so most of the interaction graph can be preserved while only a minor part has to be recomputed. Moreover, if only the weight of the rotamer library's energy term is changed, then neither the static frame energies nor the energies of pairwise interactions have to be recomputed. The parameters were optimized such that all side chains were in their crystal environment, interacting with crystal neighbors of the asymmetric unit. This is particularly important for polar side chains with contacts between asymmetric units.

Accuracy versus accessible surface area and percentile electron density

The smoothed curves in Figures 8 and 9 were calculated using kernel density estimates53 by calculating probability density estimates of correctly and incorrectly predicted side chains as a function of relative accessible surface area (RSA) and percentile of electron density:

where A_i is the RSA of residue i or its electron density percentile, and K is a Gaussian kernel with bandwidth h:

The prediction rate at A is calculated using Bayes' rule:

where p(A|corr) and p(A|incorr) are calculated using the expression for f(A) using the correctly predicted and incorrectly predicted side chains, respectively. p(corr) and p(incorr) are just the frequencies of correctly and incorrectly predicted side chains overall. Kernel density estimates exhibit unfavorable behavior at boundaries, and so the data were reflected across A = 0 to account for this. No correction was made at A = 100.

Library

SCWRL4 has been redesigned as a library, so that its optimization engine can be used in other scenarios such as protein design. It uses a delayed computation model. First, the model data (referred to as the workspace) is defined by specifying the location of amino-acid residues with appropriate rotamers via calls to functions of the library. After that the calling program uses the SCWRL4 engine in the library to derive the optimal assignment of rotamers. This will request that SCWRL4 calculate all the required energies and perform combinatorial optimization after which for each residue one of its rotamers will be marked as optimal. This information can then be used in other applications. The SCWRL4 library keeps the model alive for further usage even after the optimal assignment has been found. This means that after the optimal rotamers have been resolved, some modifications can be introduced into the model and the optimization requested again. In this case, SCWRL4 will recalculate only those energies that need to be modified due to changes in the model, while for energies between persistent objects it will use cached values.

Availability

SCWRL4 is available at http://dunbrack.fccc.edu/scwrl4.

Training and testing sets

We used separate training and testing sets of 100 and 379 proteins respectively to optimize and test various parameters and algorithmic choices for development of SCWRL4. Details of these sets are given in Supporting Information. The resolution cutoff for both sets was 1.8 Å, and the maximum mutual sequence identity was 30%. All calculations performed later are either on the asymmetric unit or using crystal symmetry, although in each case the accuracy results are compiled on sets consisting of only one chain of each sequence.

Accuracy of SCWRL4

The overall accuracy of SCWRL4 is presented in Table I for those side chains with electron density above the 25th percentile. The accuracy is reported in three different ways. First, for each side-chain type, the conditional accuracy for each dihedral degree of freedom, χ_i, is reported. This is the percent of χ_i that is correct, given that the χ angles closer to the backbone, χ_i−1,…,χ₁ are also correct. So for instance, for Met, for those residues with both χ₁ and χ₂ correct, 77.1% have χ₃ correct. The column “ALL” counts only those residue types with that degree of freedom. Second, for each side-chain the absolute accuracy at each degree of freedom is the percentage of all residues of that type such that χ_i,…,χ₁ are all correct. So for instance, for Met, 60.9% of all residues are predicted correctly for χ₁, χ₂, and χ₃. Finally, the average RMSDs are given for each side chain type. For a set of residues of the same type, there are two ways to calculate the RMSD. First, one can calculate the RMSD for each residue, and then average these values. Second, one can do the sum of square distances over all of the atoms of all of the residues, take the mean, and then the square root. The values are not the same. We use the former definition.

The numbers most frequently cited for side-chain prediction accuracy are the χ₁ and χ₁₊₂ rates over all side-chain types, where χ₁₊₂ is the absolute accuracy at the χ₂ degree of freedom in Table I. These values are 89.3 and 79.7% respectively for the 25–100th percentiles of electron density. For all side chains, these values fall to 86.1 and 74.8%. We exclude the bottom 25% because of the inherent uncertainty in these conformations (see later). For 10 of 18 residue types, the χ₁ accuracies exceed 90%, and these are predominantly the aliphatic and aromatic residue types. Ser is the most difficult to predict, with an accuracy rate of 75.8%.

In Table II, we show the improvement in prediction accuracy of SCWRL4 over SCWRL3 for each residue type and for the conditional and absolute accuracy measures. The overall improvement in χ₁ accuracy is 3.5% (SCWRL4 accuracy—SCWRL3 accuracy on the same test set of 379 proteins). The largest improvements are in Trp, Arg, Gln, Glu, Met, Asp, Asn, and Ser all of which exceed 6% improvement in average absolute accuracy.

The improvement in accuracy in SCWRL4 over SCWRL3 was achieved through a number of different changes in the sampling space, the energy function, and the algorithm. Each of these feature changes was chosen and/or optimized on the basis of improvement in the training set of 100 proteins. In Figure 6, we show the effects of each change added to SCWRL3 (boxes R through T) or each change removed from the final SCWRL4 protocol (boxes A through I). The figure demonstrates that the effect of each feature is context-dependent; that is, the effect is different when added to SCWRL3, which contains none of the new features versus when it is removed from the final SCWRL4, which contains all of the new features. The directed graph leading from SCWRL3 to SCWRL4 along the outside of the figure shows the improvements as each feature is added consecutively. The most important changes include the FRM, the new rotamer library, the addition of a hydrogen bonding function, changes in the atomic radii used, and using a larger percentage of the rotamer library (the top 98% of rotamer density, instead of 90% as used in SCWRL3). The RRM (box C) is 2.01% less accurate in average absolute accuracy than the full FRM (box A). The decrease in χ₁ and χ₁₊₂ accuracies are 1.4 and 2.8%, respectively.

Prediction of side-chain conformation in crystals

We enabled consideration of crystal symmetry in side-chain conformation prediction in SCWRL4. This is accomplished by determining the interaction graph in the context of neighboring chains to the asymmetric unit within the crystal. Thus, a residue in the graph may have a neighbor on the other side of the protein, if that residue makes contact with that residue in a crystal neighbor. Crystals were built and neighbors determined as described in previous work.44 It is of interest to determine the effect on prediction accuracy when crystal symmetry is taken into account. It should be noted that this is a bona fide prediction within the crystal, since the side chains in the neighboring asymmetric units have the same conformations as the asymmetric unit whose structure is being predicted.

In Figure 7, we show the improvement in accuracy for all side chains and for those in crystal contacts when the crystal symmetry feature is enabled. The accuracy values shown are average absolute accuracy, which are averages of the χ₁, χ₂, χ₃, and χ₄ absolute accuracies shown in Table I. Improvement occurs for all side-chain types. Among all side chains, not just those in crystal contacts, the effect is strongest for those most likely to be on the surface, in particular the longer side chains, Arg, Lys, Glu, and Gln. However, when other side chain types are in crystal interfaces, their accuracy is also strongly affected by the presence of the crystal neighbors. This is especially true for Trp and Met. The χ₁ and χ₁₊₂ accuracy in the crystal for side chains with electron densities in the 25–100% percentiles are 90.9 and 82.6%, respectively. For all side chains, the values are 87.4 and 77.1%

Prediction of side-chain conformation versus accessible surface area and electron density

Exposed side chains have fewer steric constraints and are more difficult to predict accurately. We have calculated the accuracy of predictions (within 40°) as a function of the relative surface accessibility (RSA) of side chains calculated with the program naccess.54 Both the predictions and the surface area calculations were performed in the crystal environment, so that side chains in protein interfaces in the crystal are considered buried. The results are shown in Figure 8. The accuracy versus surface area was calculated using kernel density estimates as described in the Methods section. In Figure 8, the probability density of each residue type as a function of accessible surface area is shown in magenta (multiplied by 20). These curves show that in the crystal, most residue types are predominantly buried (RSA < 40%) with maxima at 0% RSA. The only exceptions are Lys, Arg, and Glu with density modes at 30–40% exposure.

The frequency of accurate predictions is shown for each side-chain type for all side chains with RSA > 0%:χ₁ (black), χ₁₊₂ (red), χ₁₊₂₊₃ (orange), and χ_1+2+3+4 (blue). As expected, less accessible side chains are predicted more accurately than accessible side chains. Note that at high RSA, the estimates can be quite noisy due to very small counts, especially for Cys and Trp. Separate points are plotted for those side chains with 0% RSA (using the same color scheme), calculated separately from the kernel density estimates shown in the curves, For completely buried side chains in the crystal, 96.0% of χ₁ and 91.5% of χ₁₊₂ dihedrals are correctly predicted. All side-chain types are predicted with greater than 95% accuracy for χ₁ except Cys (94.4%), Pro (91.9%), and Ser (79.1%).

We have previously shown that rotameric side chains have higher electron density than nonrotameric side chains, and that nonrotameric side chains are likely to be disordered in a manner similar to side chains that are annotated in PDB files as existing in more than one χ₁ rotamer in the crystal.50 Since SCWRL4 predicts only one conformation per side chain, it seems likely that side chains with lower electron density should be harder to predict, since they may be placed in only a portion of the observable electron density. In Figure 9, we show prediction accuracy as a function of the percentile of electron density of the entire side chain, using the same color scheme as that in Figure 8. The curves are also calculated with kernel density estimates. Low percentiles correspond to low-electron density, disordered side chains, and high percentiles correspond to high-electron density, well-ordered side chains. Some low-density side chains may be incorrectly placed in the density.

For all side-chain types and all degrees of freedom, the accuracy rises with increasing electron density. For some degrees of freedom, this is only true at low electron density, and the accuracy plateaus above the 30th percentile. But for the longer side chains, the increase in accuracy extends from 0 to 100%. For all side chains, the χ₁ accuracy increases from 69.0 to 94.6% at 0 and 100th percentiles of electron density. For χ₁₊₂, the equivalent numbers are 52.4 and 82.3%. At the highest electron densities, the χ₁₊₂₊₃ and χ_1+2+3+4 accuracies are 71.9 and 56.5%.

CPU time

In Table III, we show a comparison of the CPU time for the testing set of 379 proteins for SCWRL3 and SCWRL4, for both the RRM and the FRM and for the asymmetric units and crystals. The mean, median, and maximum CPU times over the testing set reveal that the different calculations have different properties. SCWRL3 is very fast on most proteins, but on a small number of structures takes exceedingly long times. On two ASUs, SCWRL3 failed to finish and on one protein took 1409 s. SCWRL4 with the RRM model takes slightly longer than SCWRL3 with median times of 1.51 and 1.27 s, respectively. This is due to slightly longer times required for the more complicated energy function and denser and larger graphs that result. However, the maximum time for the SCWRL4 RRM is 72 s and the mean time only 4 s.

For the FRM model, SCWRL4 takes a median of 7.9 s, a mean of 12.2 s, and a maximum of 98 s. Thus, the median time is about 6.3 times slower for the SCWRL4 FRM calculation compared with SCWRL3, and the mean time is 1.5 times slower. SCWRL4 is also able to calculate the conformations of proteins in the crystal environment, taking account of crystal symmetry. Calculation with crystal symmetry takes 1.7–1.8 times that of the ASU for the FRM model of SCWRL4. The effect on the RRM model is somewhat larger.

The calculations were performed on a machine with an AMD Athlon 64 × 2 Dual Core Processor 4400+ at 2.21GHz, with 3.25 GB of RAM, and running by 32-bit Microsoft Windows XP Professional Service Pack 3.