GEM: A Gaussian evolutionary method for predicting protein side-chain conformations

Overall accuracy of prediction and comparison with other methods

The overall accuracy of GEM in predicting the side-chain conformation of 38 test proteins (4313 residues) is shown in Table 1. The geometric parameters evaluated were those used commonly by others, namely the χ₁ and χ_{1 + 2} angles for all residues and the root mean square deviation (RMSD) error in side-chain heavy atoms. The 38 structures selected constitute a minimal set encompassing most of the common structures tested in different prediction methods, allowing comparison with these other methods (Holm and Sander 1992; Dunbrack and Karplus 1993; Laughton 1994; Hwang and Liao 1995; Samudrala and Moult 1998; Looger and Hellinga 2001). The results of the comparison are shown in Tables 2 and 3.

As shown in Table 1, GEM yielded values for the χ₁ angles that were within 30° of those determined from the crystal structure in 80% of cases and side-chain atomic positions with a mean RMSD of 1.36 Å. For core residues with <20% solvent exposure, the prediction accuracy increased to 93% for χ₁ and the mean RMSD was reduced to 0.92 Å. Thus, even though no energy minimization was applied to the final structure and the bond lengths and bond angles used were those of standard templates and not of the individual residues in the individual X-ray structures (see Materials and Methods), our results for core residues, on a limited test set, were approaching those from a recent study in which a very detailed rotamer library (7560 rotamers) and experimental bonds and angles for every residue were used (Xiang and Honig 2001).

In terms of residue types, it appears easier to predict the orientation of small nonpolar or aromatic side chains, as evidenced by a significantly better accuracy in χ₁ angles (90%) being achieved for these residues (Fig. 2a,b ▶). This means that prediction errors are contributed mainly by polar and charged amino acids, for which interaction with the solvent plays a significant role in determining their side-chain orientation. The largest error was the 58% χ₁ accuracy for serine, which was probably attributable to its small size and its conformation therefore being dictated by hydrogen bonding (Koehl and Delarue 1994). These results indicate that, in general, it is harder to predict the conformation of those amino acids in which the side-chain conformation is more susceptible not only to steric hindrance, but also to environmental factors, such as hydrogen bonds, salt bridges, and solvent interactions. This observation is consistent with the findings of previous studies (Table 3; Dunbrack and Karplus 1993^; Koehl and Delarue 1994), as well as with the fact that the conformation of residues in the protein core can be more accurately predicted (Levitt et al. 1997).

In general, it is neither straightforward nor completely fair to compare the results of different protein side-chain prediction methods, as different accuracy measures, energy functions, and proteins have been used, and, with the exception of the studies of Bower et al. (1997) and Xiang and Honig (2001), tests have been performed on a rather small set of proteins (Tables 2, 3). Despite this, a number of common characteristics emerge. For example, whereas no obvious structures or classes of proteins emerge as particularly easy or difficult to predict by any of the different prediction methods, on the whole, the prediction accuracy does not differ significantly from one method to another. Furthermore, as mentioned above, the prediction accuracy for χ₁ angles in terms of amino acid type is quite consistent between these methods (Table 3). The availability of SCWRL (Side-chain placement With a Rotamer Library) (Bower et al. 1997) through the Internet (http://www.fccc.edu/research/labs/dunbrack/scwrl/), one of the more recent and widely used side-chain modeling programs allowed a more straightforward comparison, and the results of SCWRL we obtained on the 38 test proteins using the same accuracy criteria of GEM (Table 1) were also similar, with the overall average RMSD being 1.46 Å, χ₁ accuracy being 81%, and χ₁₊₂ accuracy being 64%. These results suggest that the accuracy of GEM is comparable with those of previous prediction methods. However, the GEM approach, as discussed below, can be used to analyze elements of methodology, such as search scheme and energy function, and should therefore help in moving toward error-free prediction of protein side-chain conformations.

Evaluation of the energy function used

The main objective of this study was to evaluate whether the evolutionary algorithm that we recently developed and applied to flexible ligand docking (Yang and Kao 2000a) was also applicable to the prediction of side-chain conformations of proteins. To simplify the task, we adopted a typical van der Waals–type energy function used in previous studies (see Materials and Methods). However, to enable the global continuous search mechanism of GEM to energetically discriminate between different torsions, we found it necessary to add a torsion term. In two previous studies in which a very large set of rotamer states were sampled, a torsion function was also used (Lee and Subbiah 1991; Xiang and Honig 2001). The torsion term used in the present work was that employed in the HIV protease-inhibitor docking study of Gehlhaar et al. (1995).

The fact that although we did not attempt to refine any parameters of the energy function used, we still achieved a comparable prediction accuracy attested to the viability of the use of GEM for protein side-chain prediction. However, with uncertainty in the scoring (fitness) function, the robustness of GEM was difficult to assess. To address this question, we made use of the high adaptability of GEM and simply replaced the force-field energy function with a perfect-scoring function (i.e., one that would produce zero RMSD in atom positions). As shown in Table 1 and Figure 2, c and d ▶, using the RMSD scoring function (equation 1; Materials and Methods), GEM could indeed approach perfect prediction. Not only were almost all χ₁ angles predicted to within 30° (Fig. 2c ▶; Table 1), but the absolute values of the angles were very accurately reproduced, with >90% having an error of <5° compared to a value of only 30% when the force-field energy function (equation 2) was used (Fig. 3a ▶). The residual errors can be attributed to the use of not completely flexible side-chain templates (Materials and Methods), which have rigid bond lengths and bond angles. This attribution is supported by the work of Xiang and Honig (2001), which showed that the use of a standardized geometry could lead to considerable errors. This was most evident in the RMSD results for lysine and arginine, and, to a lesser extent, those for glutamate and glutamine (Fig. 2d ▶), that is, those amino acids with side chains that are charged or polar and are relatively long and flexible. Interestingly, in the case of lysine and arginine, although almost all the χ₁ angles were predicted without error, the RMSD value was as large as that obtained using the force-field energy function (Fig. 2 ▶), reinforcing the fact that either RMSD or χ₁ alone cannot provide a complete assessment of side-chain prediction accuracy. It is also worthy of note that GEM converges much faster with the perfect-fitness function (Fig. 3b ▶).

Limitation of discrete search

The purpose of discrete search, such as the use of a rotamer library or uniformly divided torsion angles, is to reduce the side-chain conformational space to render computational search tractable. However, it frequently happens that many side-chain conformations cannot be covered if only discrete angle values are used. Indeed, the average prediction accuracy of GEM with discrete search alone was 70.9% for χ₁ and 1.72 Å for RMSD for the 38 test proteins (Table 4). GEM with both discrete and continuous search consequently enjoyed an improvement of 9.1% in average χ₁ prediction and 0.36 Å in average RMSD of side-chain positions (Fig.2a,b ▶). Furthermore, even using the perfect-fitness function, GEM with discrete search alone resulted in substantial errors, especially for flexible residues, such as serine, lysine, and arginine (Fig. 2c,d ▶).

The shortcoming of discrete search over a limited number of rotamer states (103 in this study) is even more evident in the case of strained conformations (Schrauber et al. 1993), which may be defined as those with a χ₁ angle deviating by >30° from the three energetically favored values (+/−60° and 180°). As shown in Figure 4 ▶, the discrete search with GEM yielded a difference as large as 30% in the χ₁ accuracy between nonstrained and strained side-chain conformations. Such a large difference could often be easily overlooked using statistical averages because strained residues represent a small minority of all amino acid side-chain conformations (e.g., of the 4313 residues examined here, only 177 [5%] were strained).

A more detailed rotamer library may improve side-chain prediction accuracy (De Maeyer et al. 1997; Liang and Grishin 2002). Using a very large library with 7560 rotamer states, Xiang and Honig (2001) indeed obtained very good performance, especially for core residues. However, larger rotamer states have not generally proven more accurate for side-chain prediction than smaller states (Holm and Sander 1992; Laughton 1994; Väsquez 1995). Using a library of 330 states (De Maeyer et al. 1997), GEM with the energy-based scoring function (equation 2) yielded very marginal improvements over the results with 103 states when only the discrete-search mechanism was used, and essentially no improvements when discrete search was combined with continuous search (Table 4). However, it should be noted that because the statistics of this larger library were not available, the statistics of the 103 states were transferred and equal probability was assigned to each of the additional states in the larger library, and, as such, the calculation could not take full advantage of the evolutionary parameters optimized in GEM. Interestingly, using the exact RMSD scoring function (equation 1), the discrete-only prediction of GEM improved significantly with the larger library (χ₁ from 78.8% to 85.6% and RMSD from 1.26 Å to 0.87 Å; Table 4). These results may suggest that the scoring function is a key accuracy-determining factor in side-chain prediction, and that the efficiency of utilizing a larger library depends on search mechanisms. The suggestion is consistent with a very recent finding of Liang and Grishin (2002), who specifically optimized a scoring function to achieve very accurate side-chain predictions.

In summary, we have demonstrated the robustness and adaptability of GEM for exploring the conformational space of protein side chains and efficiently finding the combinatorial solution under the constraint of the fitness function used. The key novelty of the present work is the seamless ability of GEM to blend global discrete search and global continuous search, the former being required for efficiency and the latter for accuracy, and to allow them to work cooperatively. This was achieved through the incorporation of a number of genetic operators, each having unique search mechanisms (Table 5). For example, whereas the rotamer mutation operator performs a discrete search on the global rotamer space of the rotamer library, the self-adaptive Gaussian mutation performs a local, but continuous, search on torsional conformations, etc. Importantly, the flexibility of GEM should allow us to begin to systematically improve the forms and parameters of energy function for protein side-chain and other protein-structure prediction problems.

GEM parameters and computational details

The GEM parameters used in this paper are listed in Table 6. These parameters were selected after many attempts to predict conformations for test proteins with various initial values. GEM optimization stops when either the convergence is below a certain threshold value or when the iterations exceed a preset maximum value. In this paper, the maximal number of generations was set to be 100 + K/2 or 250 if 100 + K/2 > 250, where K is the residue number of a protein. A set of 38 high-resolution crystal structures of proteins (resolution better than 2 Å) ranging in size from 46 to 325 amino acid residues was used to test the performance of GEM. This set was selected from those used in many previous studies (Holm and Sander 1992; Dunbrack and Karplus 1993; Laughton 1994; Hwang and Liao 1995; Samudrala and Moult 1998; Looger and Hellinga 2001) to compare our results with those obtained using other methods. All calculations were performed on a 500-mHz Pentium III processor. A typical run time for a protein of 300 amino acids was ∼25 min.

Root mean square deviation (RMSD) of atomic positions and the percentage of side-chain dihedral angles that were correctly predicted within 30° were used to assess the accuracy of the prediction. The RMSD was calculated using the formula

(1)

where M is the number of atoms in a protein and (X_i, Y_i, Z_i) and (x_i, y_i, z_i) the coordinates of the ith atom of the X-ray crystal and the predicted structure, respectively. The calculation of RMSD included the heavy atoms of side chains. The accuracies for χ₁ and χ₁₊₂ were calculated excluding glycine and alanine residues. The accuracy for χ₁₊₂ refers to those cases in which both χ₁ and χ₂ were correctly predicted.

Rotamer library and side-chain construction

A main-chain independent rotamer library was built using a modification of the method of Tuffery et al (1991). For each amino acid, only up to ten of the highest populated rotamers were considered. The number of rotamers was 1 for Pro; 3 for Cys, Ser, Thr, Val, Asp, and Phe; 4 for Asn; 5 for Ile; 6 for His, Leu, and Tyr; 7 for Trp; and 10 for the remaining amino acids (excluding Ala and Gly, which have no rotamers). The total number of rotamers was 103.

The side-chain atoms were geometrically constructed by placing side-chain templates onto the main chain of the X-ray structures, using the side-chain dihedral angles generated by GEM. These templates have standard bond lengths and bond angles according to the AMBER force field (Weiner et al. 1984). The initial dihedral angles of the side chains came either from the rotamer library or from an assigned value within the feasible region (−π to π). GEM then adapted these dihedral angles to search for optimal side-chain conformations by minimizing the scoring function.

Scoring energy function

Our scoring energy function was modified from Levitt (1983) and Hwang and Liao (1995):

(2)

where E_vdw is the van der Waals interaction potential, E_tor the torsional potential, and E_Hbond and E_Sbond are the potentials of hydrogen bonds and disulfide bonds, respectively. E_vdw was computed using a Lennard-Jones 6–12 potential:

(3)

where ɛ_ij and r_ij are constants which depend on the chemical characteristics of atoms i and j; R_ij is the distance between the atoms i and j; M is the number of atoms in a protein; NB is the number of atoms within a preset distance (8 Å) of atom i; ɛ_ij is the depth of the energy well, and r_ij is the equilibrium interatomic distance for the van der Waals interaction between atoms i and j. The same function form was used to compute the potential of the hydrogen bonds (E_Hbond) and disulfide bonds (E_Sbond). Table 7 shows the values used for these parameters. We restricted the E_vdw energy to a maximal value of 20 kcal/mole for each atom pair, as described by Levitt (1983), to avoid infinite energies.

For the torsional energy, the equation of Gehlhaar (1995) was used:

(4)

where K is the residue number of a protein and chi is the number of dihedral angles of a residue. The values for A, n, and τ₀ were 3, 3, and π, respectively, for the sp3-sp3 type and 1.5, 6, and 0, respectively, for the sp3-sp2 type (Gehlhaar et al. 1995).

GEM algorithm details

Here, we provide an outline of our evolutionary approach for predicting protein side-chain conformations, which can be represented by adjustable variables of dihedral angles as

(5)

where n is the number of dihedral angles of a side-chain conformation. Generally, the steps involved are as follows:

1. Initialize the side-chain conformation of each residue on a given backbone. The initial values for the dihedral angles are selected either from the rotamer library or from the feasible region (−π, π). Repeat this N times to generate the initial population of N side-chain conformations for a protein structure. Evaluate the objective value of each conformation based on the scoring function.

2. Change the dihedral angles of side-chain conformations by genetic operators to generate offspring. Evaluate the objective values of the offspring.

3. Use selection operators to select N solutions from the side-chain conformations of both parent and offspring solutions.

4. Repeat steps 2 and 3 until one of the terminating conditions is satisfied.

Main procedure

In the following subsections, we present the details of our approach for side-chain prediction. The method integrates a global discrete-search mechanism based on a rotamer library and continuous-search mechanisms based on Gaussian mutations. The basic structure of the method is as follows (Fig. 5 ▶): N solutions are generated as the initial population. Each solution (side-chain conformation) is represented as a set of three n-dimensional vectors (θⁱ, σⁱ, υⁱ), where n is the number of adjustable variables (dihedral angles) of a side-chain conformation and i = 1,. . .,N. The vector θ in equation 5 represents the adjustable variable to be optimized, and σ and υ are the step-size vectors of decreasing-based mutation and self-adaptive Gaussian mutation, respectively (Yang and Kao 2000a,b). In other words, each solution, θ, is associated with some parameters for step-size control. In this paper, the initial value of θ_jwas randomly selected either from the rotamer library or from −π to π in radians. The initial step sizes, σ and υ, were 0.8 and 0.2 radians, respectively.

The main optimization procedure consists of three stages in one generation: a rotamer-search stage, a decreasing-based Gaussian mutation stage, and a self-adaptive Gaussian mutation stage. The rotamer-search stage is a discrete-search mechanism that uses the rotamer mutation operator to find an optimal combination of rotamer conformations. The decreasing-based Gaussian mutation and self-adaptive Gaussian mutation are continuous-search mechanisms that mutate dihedral angles to find an optimal side-chain conformation in continuous-search spaces. As shown in Figure 5 ▶, each stage uses a general procedure, FC_adaptive, with two parameters to generate a new quasi-population (with N solutions) as the parent of the next stage. These stages differ only in the mutations used. The recombination and mutation operators will be described below.

The FC_adaptive procedure employs two parameters, the population operator (P, with N solutions) and the mutation operator (M), to generate a new quasi-population. The main purpose of the FC_adaptive procedure is to produce offspring and then perform the family competition. Each individual in the population sequentially becomes the "family father." This family father and another solution randomly chosen from the rest of the parent population are used as parents for a recombination operation, with a probability of recombination of p_c; then a mutation operates on the new offspring or the family father (if recombination does not occur). For each family father, this procedure is repeated L times (family competition length). Finally, L children are produced, but only the one with the lowest objective value survives. Because we create L children from one family father and perform a selection, this is a family-competition strategy.

For easy description of operators, we use a = ( θ^a, σ^a, υ^a) to represent the family father and b = ( θ^b, σ^b, υ^b) as another parent (only for the recombination operator). The offspring of each operation is represented as c = ( θ^c, σ^c, υ^c). We also use the symbol θ_j^d to denote the jth dihedral angle of a side-chain conformation d.

Recombination operators

Modified discrete recombination

Because experience indicated that our method was more robust if the child inherited genes from the family father with a higher probability, we therefore modified the original discrete recombination (Bäck 1996) as follows:

(6)

Intermediate recombination

This operator is applied in the continuous-search stages as follows:

(7)

(8)

where ω is σ or υ, depending on the mutation operator applied. For example, if the self-adaptive Gaussian mutation was used in this FC_adaptive procedure, ω is υ. β is a constant set as 0.5 in this work.

Mutation operators

Mutations are the main operators of our method. After recombination, a mutation operator is applied to the family father or to the new offspring generated by a recombination operator.

Rotamer mutation operator

This operator is used at the rotamer-search stage to find a combination of rotamer conformations. For each residue, this operator is biased to select rotamers of higher probabilities and mutates all of the dihedral angles of a residue according to the rotamer library. For example, this operator changes 3 dihedral angles (θ_j, θ_{j + 1}, and θ_{j + 2}) if the residue is Gln, Glu, or, Met. For each of these dihedral angles, this operator is applied with probability 0.2 as follows:

(9)

where γ_ki and p_ki are the angle value and probability, respectively, of the ith rotamer of the kth residue type. The values of γ_ki and p_ki are defined in the rotamer library.

Self-adaptive Gaussian mutation

The mutation is accomplished by first mutating the step size υ_j and then mutating the dihedral angle θ_j as follows:

(10)

(11)

where N(0, 1) is the standard normal distribution and N_j(0, 1) is a new value with distribution N(0, 1) that must be regenerated for each index j. We followed Bäck (1996) in setting τ and τ as (√2n)⁻¹ and (√2√n)⁻¹, respectively.

Decreasing-based Gaussian mutations

The decreasing-based Gaussian mutation is accomplished by mutating the step-size vector σ with a fixed decreasing rate γ = 0.95 as follows:

(12)

(13)

Previous results (Yang and Kao 2000b) showed that self-adaptive mutations converge faster than decreasing-based mutations, whereas, for rugged functions, self-adaptive mutations tend to yield optimization results that are confined to local minima more easily than decreasing-based mutations.