3264
J. Phys. Chem. 1996, 100, 3264-3272
Minimizing Reduced-Model Proteins Using a Generalized Hierarchical Table-Lookup Potential Function John R. Gunn De´ partement de Chimie, UniVersite´ de Montre´ al, C.P. 6128, Succ. Centre-Ville, Montre´ al, Que´ bec H3C 3J7, Canada, and Centre de Recherche en Calcul Applique´ , 5160 boul. De´ carie, Bureau 400, Montre´ al, Que´ bec H3X 2H9, Canada ReceiVed: August 21, 1995; In Final Form: NoVember 7, 1995X
A generalized multidimensional contact potential is introduced for simulating a simplified model of protein tertiary structure. The contact energy between two residues is a simultaneous function of five interatomic distances, which allows orientational correlations between the backbone and side chains to be represented. The minimization of this potential is enhanced by the use of three levels of coarse-grained potentials as screening functions which can cheaply reject unfavorable structures. The approximate potentials are systematically derived from the full contact potential by averaging over the unwanted degrees of freedom. Results are shown for the minimization of three medium-sized (60-146-residue) proteins and compared with earlier results. The performance of the hierarchical algorithm as a function of the coarse-graining scheme and screening cutoffs is discussed, as well as the dependence of the potential evaluation on system size. The results from the three test proteins show a significant improvement from the N2 scaling observed without the use of the screening potentials.
I. Introduction The ab initio prediction of protein tertiary structure continues to be one of the outstanding problems in computational chemistry. From the point of view of computer simulation, the problem is one of minimizing a generally complicated 3Ndimensional function, namely, the effective free energy as a function of the coordinates. This presents two formidable hurdles: developing a sufficiently accurate representation of the energy function and then effectively sampling it in order to find the global minimum.1 In this paper we present a general algorithm which makes use of a recently developed hierarchical structure2-4 to significantly improve the sampling efficiency of a Monte Carlo simulated annealing calculation with an increased specificity of the potential function. Some implications for the development of improved potential functions will also be discussed. The basis of the hierarchical algorithm consists of a multilevel coarse graining of the potential function which allows the trial conformations to be evaluated with varying degrees of resolution. As a computational technique, the idea of hierarchical resolution is well-known in the context of grid-based equation solvers such as those found in computational fluid dynamics.5 The underlying goal is to reduce the size scaling of the computational effort by shifting most of the work to a faster low-resolution representation while still maintaining the final accuracy of the high-resolution representation. In the case of protein structure, the simplest coarse graining of the chemical structure is at the level of the primary sequence, namely, the constituent amino acid residues, and this is the starting point of the current model. A further level of coarse graining arises naturally at the level of the secondary structure, R-helices and β-strands, where there is to some extent a separation of energy scales between the specific local interactions which stabilize the secondary structure and the nonlocal packing forces which assemble the structural units into the folded structure.6,7 We therefore address the problem of X
Abstract published in AdVance ACS Abstracts, January 15, 1996.
0022-3654/96/20100-3264$12.00/0
determining tertiary structure with a specified secondary structure and thus only consider a potential function which represents the latter type of interaction which we further assume can be reasonably approximated as a potential between the distinct structural units. This is the basis of the “cylinder and sphere” model2 in which the secondary structural units are treated as rigid bodies subject to a coarse-grained potential.8 The combination of the residue-based and cylinder-sphere models leads to a truly hierarchical algorithm in which the trial moves in the Monte Carlo minimization of the residue-based potential are themselves the result of a minimization of the cylinder-sphere potential. Although the hierarchical algorithm has proven effective at sampling the gross features of folding topology, it has certain deficiencies which will be addressed here. One fundamental difficulty arises from the use of distinct physical models at each level of resolution. The cylinder-sphere model is inherently limited in the accuracy with which it can approximate the residue-based potential, in particular the excluded volume which has the form of a solid core. This means that there is no way to systematically refine the coarse graining in order to improve the convergence of the minimization. This also applies to other terms in the potential energy such as contact potentials for which there is no suitable representation in the cylinder-sphere coordinates. For this reason, a refinement algorithm was developed4 in which all trial moves are carried out with the same representation of the molecule. This also removes the restriction of having the trial moves consist of substitutions of entire loops, since the atomic coordinates are moved directly at each step. The greatest limitation, however, is in the simple form of the potential function. Although the hydrophobicity potential9 was effective at generating compact structures, it was insufficient to distinguish the folded topology of the native from other plausible misfolded structures. In order to develop a more accurate energy calculation, a generalized contact potential in the form of a lookup table was introduced.3 The added complexity, however, means that in a straightforward minimiza© 1996 American Chemical Society
Minimizing Reduced-Model Proteins tion of the generalized potential the computational effort is completely dominated by the evaluation of the energy. The motivation of the present work is therefore to demonstrate that a hierarchical coarse-graining strategy can be used to effectively minimize such a generalized tabular potential with a substantial improvement in computational efficiency. In this way, the numerical advantage which motivated the original cylindersphere model can be partially recovered while removing some of the limitations on the final resolution of the model. II. Simulation Methodology The conformation of the protein molecule is represented by a series of conformational states, one per residue, which specify the values of all internal degrees of freedom of a residue. In the present model, the only variables are the two backbone dihedral angles, φ and ψ, with all other coordinates held constant for all states. There is a discrete set of possible states, chosen to span the accessible region of the Ramachandran map and which are the same for all amino acid types. The states lie on alternate sites of a five degree square grid (648 points) in φ-ψ space, restricted to the populated regions described by Rooman et al.10 leaving 532 points. The geometric model of the structure is thus sequence-independent, and includes only those atoms determined by the backbone conformation (with the β-carbon being fictitious in the case of glycine), and can be completely specified by a string of integers specifying the states of each residue. The trial moves in the minimization consist of changing the states of one or more residues, resulting in a rigid-body rotation of the remainder of the molecule. The secondary structure is imposed simply by initializing specified residues to states corresponding to typical helix or strand geometries and holding them fixed throughout the calculation. The potential function has the form of a generalized contact potential between residues. In principle, six coordinates are necessary in order to completely specify the relative orientations in space of two residues, and thus six interatomic distances can be used to effectively take into account all interactions between rigid structures. Alternatively, the four distances between Rand β-carbons on different residues can be thought of as a cylindrically symmetric approximation with the axes corresponding to the backbone-side chain direction.3 An intermediate value of five was chosen for the present work, namely, |rCβ,i - rCβ,j|, |rCβ,i - rC′,j|, |rCβ,i - rN,j|, |rC′,i - rCβ,j|, and |rN,i rCβ,j|. This was due to storage constraints and also to match the form of the empirical potential, as described below. The N, C′, and Cβ atoms are used as reference points, with the CR atoms being necessarily located at the pyramidal positions in the middle. Besides the Cβ-Cβ distance, the potential makes use of the distances between the Cβ atoms and the other backbone atoms to provide an indication of the backbone orientation. As mentioned above, one other distance would be needed to uniquely determine all atomic positions, and any more than that would be redundant. The potential is thus a fivedimensional lookup table, with in principle an independent value for each relative orientation of each amino acid pair, with a resolution determined by the discretization of the distance indices. The potential function consists of two principal components, an empirical packing potential and some form of excluded volume. The empirical potential used here is a pairwise contact potential developed by Bryant and Lawrence.11 This potential has the general form of a distance-based table with interactions between side chain and backbone “centers” which are defined as points 2.4 Å from the CR atoms in the Cβ direction and as points midway between adjacent CR atoms, respectively. This
J. Phys. Chem., Vol. 100, No. 8, 1996 3265 potential was chosen because the two types of interactions incorporate some information about the relative orientations of the side chains and backbone since, although the potential itself is pairwise additive, the parameters in the table were optimized simultaneously. In general, however, the tabular potential is not separable and is a function of all coordinates. This is especially true for the excluded volume where the positions of all other atoms are extrapolated from the relative orientations of the residues. This is a general property of the coarse-graining scheme, namely, that a renormalization of pairwise interactions leads to many-body correlations in the remaining degrees of freedom. The excluded volume plays an important role in the minimization because there is nothing in the contact potential to prevent the molecule from collapsing. This is because the parametrization was based on known structures in the Protein Data Bank which automatically have a reasonable density, and therefore the potential lacks a realistic density dependence. On the other hand, a steeply repulsive hard core implies enormous barriers in the potential surface which hinders the minimization. It is therefore necessary to treat the hard-core repulsion as a separate constraint which is annealed in during the course of the simulation. This is accomplished by assigning inaccessible relative conformations a positive energetic penalty which is an adjustable parameter in the simulation. The entries in the lookup table corresponding to “overlapping” conformations are then updated during the minimization. An additional problem arises due to the fact that the lookup table is limited to interactions within a cutoff distance and is thus slow to generate compact structures since in a randomly generated initial configuration there is little energetic gain initially in trying to fold. To alleviate this, an additional term is added to the potential directly proportional to the radius of gyration and is gradually removed during the simulation. The average density of the final ensemble is thus somewhat dependent on the balance between the overlap penalty and this “external pressure”. The minimization is carried out using a standard Monte Carlo simulated annealing algorithm, with the temperature decreasing as an exponential function of the iteration count. In addition, the overlap penalty and pressure term are varied linearly with the interaction count. The minimization is carried out simultaneously on an ensemble of 128 independent conformations. The trial moves consist of the replacement of three consecutive residue conformations, where the trial segments are selected from a precalculated list. The list is generated by randomly generating segment conformations and then screening them according to the net relative rotation of the end points, which is discretized into bins in the internal coordinate space. Only those segments falling into bins which are populated in the existing ensemble of structures are kept. The list is then sorted so that the trial moves can be selected from the same rotational bin as the segment to be replaced. This establishes a limit on the global conformational change introduced by the new segment, and thus increases the acceptance rate, while still having large enough bins to allow efficient sampling of the conformational space. In this work, a bin size of 45° was used for a total of 8192 bins in the five-dimensional conformational space. The number of trial segments in the list was 16 times the number of trial moves per molecule, or equivalently oneeighth the total number of trial moves, with typically on the order of 10 in each occupied bin at any given time. The presorting of the segment list has a minimal computational cost since the same list can be used for all molecules in the ensemble, and furthermore the time required to generate the list is small
3266 J. Phys. Chem., Vol. 100, No. 8, 1996 compared to the evaluation of the potential. The details of this procedure can be found in ref 4. III. The Hierarchical Potential The lowest-level potential function consists of a five-index lookup table as described above. The distances were discretized into 15 bins, according to the following scheme: one bin at 0-3 Å, four 0.5 Å bins from 3 to 5 Å, and ten 1 Å bins from 5 to 15 Å. This is similar to the scheme used in ref 11 which consisted of one bin at 0-5 Å and five 1 Å bins from 5 to 10 Å, but the correspondence is not exact for two reasons. First, the points used in the ref 11 potential do not correspond to real atoms and must be extrapolated from the interatomic distances used in the five-index table. This causes the bins to overlap somewhat with the expansion to larger distances. Second, the positions of the peptide centers depend on the backbone conformation, and thus combining the interactions into a single residue-residue potential represents an average over allowable values of φ and ψ. The values in the table were thus determined by taking an average over a large number of randomly generated conformations of two residues, each represented by independent fragments of three CR atoms, i.e., acetyl- and methylaminecapped alanine (also known as alanine dipeptide). The fragments are independent since the potential only includes residue pairs separated by four or more in the sequence. In each case, the bin in the five-index table was determined from the interatomic distances, and the energy was calculated according to ref 11 with the side chain-peptide terms divided by two to prevent overcounting. In addition, distances were calculated for all pairs of atoms between the two fragments (12 atoms each, neglecting saturated hydrogens), and the entire configuration was rejected if there were any van der Waals overlaps, the average therefore only including the remaining conformations. This was continued until the number of occupied bins was roughly steady, with the assumption being that empty bins were either geometrically impossible due to the connectivity, mutually intersecting, or so close to being one or the other as to have a negligible number of structures contribute to the average. Since the residue conformations are sequence-independent, the energies were calculated simultaneously for all 210 amino acid pairs from the same distribution. This averaging procedure thus replaces the excluded-volume constraints of 144 atom-atom pairs with a correlated function of the five indices, such that the labeling of a bin as accessible or not depends on how the entire fragments “fit together” and not just the calculated distances. For the empirical potential, it creates a residue-residue potential from a sum of pairwise terms. In the present case, this would seem to be an unnecessary complication; however, it is a general procedure which could be applied to a more detailed potential function where, as in the case of the excluded volume, it represents a substantial reduction in the number of degrees of freedom. The generation of the table is thus itself the first level of coarse graining within the hierarchical framework. This averaging also has the effect of smoothing out the potential surface, although the net change in the total energy of the molecule compared to the original potential is generally small. The next level of coarse graining is the reduction of the fiveindex table to a single-index residue-residue potential, based on the CR interatomic distances. This is done by taking a weighted average over all values of the other distances of conformations falling in each bin. The average is weighted in this case to better represent the minima of the full potential within bins which cover a wide range of energies and also as a means of including the overlap penalty. Since this penalty varies
Gunn TABLE 1: Average Radius of Gyration and Overlap Energy for Calculated Ensembles 1MBO
1CTF
1R69
14.9 16.3 1 6.8 5.0
10.4 12.0 2 10.9 1.0
9.8 11.1 1 3.1 3.0
native Rg average Rg native number of overlaps average number of overlaps energy per overlap
during the simulation, the fraction of overlaps in each bin is stored along with the average of nonoverlapping energies so that the total weighted average can be recalculated. The coarsegrained energy thus has the form
E)
A + e-β B + e-β
(1)
where
A)
1
∑i ′Eie-βE
(2)
∑i ′e-βE
(3)
i
Noverlap
and
B)
1
Noverlap
i
are the two fixed parameters used to specify each coarse-grained bin. Here, is the value of the variable penalty, β is the inverse temperature used for the canonical average, and the sums are over the nonoverlapping structures in each bin. The value of β was treated as an independent parameter and chosen to give the best agreement between the full and approximate potentials and was found to be largely independent of both and the annealing temperature used in the minimization. The values of the potential for individual contacts are typically on the order of 1 unit, with the overlap energies being slightly higher (see Table 1). The value of β-1 used here was 1.0, reflecting this energy scale. Higher levels of approximation were calculated as potentials acting between multiresidue segments of the chain. The segments are defined by including only every nth CR in a distance matrix and interpolating the positions of the skipped residues. The relative orientation of two segments can be defined by the positions of their end points, leading to a potential which is a function of six coordinates, consisting of the distances |rCR,i - rCR,i+n|, |rCR,j - rCR,j+n|, |rCR,i - rCR,j|, |rCR,i+n - rCR,j|, |rCR,i - rCR,j+n|, and |rCR,i+n - rCR,j+n|. If the segments were assumed to be linear, this would be sufficient to determine all interresidue distances; however, for arbitrary segments this approximation does not justify the effort of calculating six-index tables for each pair. The simplest approximation would be a single-index potential acting between the remaining residues, where the interaction at each is a sum of all the amino acids closest to it. This is a plausible approximation to the hydrophobic energy but fails badly at estimating the number of overlaps since it would count n2 each time two CR’s were too close and zero otherwise. A compromise is used in the present work, where the potential is a two-index function acting between a point and a pair of consecutive points defining a segment. The total contribution of each residue pair is determined by assigning its energy proportionally to the interactions among its neighboring segment end points. For a residue i and two adjacent residues j and k, where k ) j + n, this averaging takes the form
Minimizing Reduced-Model Proteins l1max
E(rij,rik) )
J. Phys. Chem., Vol. 100, No. 8, 1996 3267
l2max
∑ ∑
l1)l1min l2)l2min
w1w2El1,l2(c1rij + c2rik)
(4)
where the summation goes from l1min ) i - n + 1 to l1max ) i + n - 1 and from l2min ) j - [(n + 1)/2] + 1 to l2max ) k + [(n + 1)/2] - 1. The weight factors are
w1 ) 1 -
|l1 - i| |2l2 - j - k| and w2 ) 1 n 2n
(5)
and the distance coefficients are
c1 ) 1 -
l2 - j l2 - k and c2 ) 1 + n n
(6)
The potential between the point i and the segment j - k is thus a weighted sum of interactions between all of the pairs of residues within a range of 2n on each side. The distance used to evaluate each potential is also a weighted average of the two distances which serve as the reference. This also means that each pair of residues will contribute eight times to the total energy due to the fact that the ranges of summation for the segments overlap (with the sum of the weights being unity). The total energy of a given pair of residues is therefore an average of terms evaluated at eight different distances corresponding to averages of the pairs of distances associated with each segment-segment energy. This is clearly not an accurate model of the real spatial relationships of the residues but is intended simply to smooth out the potential in such a way that the total error in the energy of a molecule is more evenly distributed. In particular, the error in counting overlaps should be reduced to O(n). The total energy of the molecule is calculated by carrying out this asymmetric approximation between every nth residue and all nonneighboring pairs which define n-residue segments and then dividing by two. It should be pointed out that while the construction of the approximate potential table is highly redundant, the subsequent evaluation of the approximate potential proceeds by an ordinary table look-up, as with the original potential, using the sparse distance matrix of every nth residue. The higher-level approximations are constructed by using the first-level approximation as the residue-residue potential, which already includes the relative contribution of the overlap energy as a function of distance. The structure of the hierarchical approximations to the potential is illustrated schematically in Figure 1. This differs from the cylinder-and-sphere model in that the segments are not well-defined and the approximate potentials are necessarily extremely crude. The correlations between the approximate potentials and the full potential are shown in Figure 2 for a representative number of structures taken from an intermediate point in a typical minimization of myoglobin as discussed in section V. It is clear that the approximation can at best only give a rough idea of the quality of a structure and that it would not be productive to minimize the approximate potential independently as was the case with the cylindercylinder potential of ref 2. The approximate potentials are therefore used only as filters which serve to screen out the trial conformations anticipated to have the lowest chance of acceptance, while passing all others to the next level of evaluation with the final acceptance made only after the full potential is evaluated. This is similar in principle to a related procedure12 which can be implemented exactly for a model where there is a short-range cutoff, and the screening potential can be used to rigorously rule out any possible interactions between distant
Figure 1. Schematic drawing of the levels of approximation showing the interatomic distances which constitute the degrees of freedom at each step. Atomic coordinates used in each approximation are shown with level 0 corresponding to the full potential and level 3 to the crudest approximation. The interatomic distances which are calculated are shown as dashed lines.
segments. It was shown that the hierarchical algorithm in that case leads to a linear scaling of computational effort with chain length. In order to accommodate long-range forces, there is necessarily a tradeoff between increased efficiency and the premature rejection of trial moves which would otherwise be accepted. The practical consequences of this will be assessed below. The performance of the hierarchical algorithm clearly depends on the length of the segments used in the approximate potentials at each level along with the number of levels. For the preliminary testing, two multiresidue levels of approximation were used corresponding to n ) 3 and n ) 9. Choosing the larger segment to be a multiple of the smaller allows the same end points to be reused. This choice corresponds to a reduction of the number of variables by roughly an order of magnitude at each level. Additional higher levels would be larger than the typical size of secondary structural units and therefore unlikely to provide a useful estimate of the energy. The implementation of the screening algorithm requires that a cutoff value be established at each level of approximation, such that trial moves with approximate energies above the cutoff can be summarily rejected. This cutoff is relative to the previous energy of the structure and depends on the typical error in the approximation and the desired safety margin. Since the correspondence between the numerical values of the different approximations is generally unknown and varies with the overlap energy during the simulation, it is determined empirically by calculating the energy at all levels of approximation for each member of the current ensemble of structures. The cutoffs are then set to the lowest values that would allow all members of the ensemble to pass. For a sufficiently large ensemble, this should be a reasonable estimate of the range of values likely to be encountered in accepted moves. An additional safety margin is then added which is a user-defined constant for each level of approximation. These values can then be adjusted to control the overall performance of the algorithm as discussed in more detail in section V. A schematic outline of the algorithm is
3268 J. Phys. Chem., Vol. 100, No. 8, 1996
Gunn
Figure 3. Schematic outline of the implementation of the hierarchical potential, showing the construction of the approximate tables and the screening of the trial moves.
Figure 2. Correlations between full and approximate potentials, shown as a contour map of the distribution of 655 360 sample trial moves. Contours are at relative densities of 0.01, 0.1, and 0.5. (a) Singleindex residue-residue potential; (b) three-residue segment-segment potential; (c) nine-residue segment-segment potential. The correlation coefficients of least-squares-fit lines are respectively (a) 0.96, (b) 0.78, and (c) 0.50.
presented in Figure 3, showing both the initialization and evaluation stages of the hierarchical potential. A typical example of the variation in the empirical energy shifts during a minimization is shown in Figure 4. The number of structures evaluated at each level of approximation is also shown for the same run in Figure 5. Each iteration consists of 1024 trial moves for each of 128 structures in the ensemble. Even though there is significant fluctuation in the calculated cutoffs, the overall screening is relatively stable with time. Figure 5 shows a short initial equilibration period followed by a longer decay in the number of full potential evaluations as the acceptance
Figure 4. Energy shifts applied to the potential approximations as a function of iteration count: (a) level 3 approximation (solid line), (b) level 2 approximation (long-dash line), and (c) level 1 approximation (short-dash line).
rate decreases with the annealing temperature. In this example, the acceptance rate of the structures evaluated with the full potential fell smoothly from 70% to 10% during the simulation. IV. Results The algorithm was tested on the following three examples of small globular proteins: sperm whale myoglobin (1MBO13), 146 residues, 8 helices, and 25% loop; L7/L12 ribosomal protein C-terminal domain (1CTF14), 66 residues, 3 helices, 3 strands, and 26% loop; and phage 434 repressor N-terminal domain (1R6915), 60 residues, 5 helices, and 33% loop. The value of the overlap energy was adjusted in order to generate an ensemble of structures with a reasonable average radius of gyration at the end of the simulation. The overlap energy is turned on linearly, with the initial and final values differing typically by
Minimizing Reduced-Model Proteins
Figure 5. Number of structures evaluated at each level of approximation for the same simulation as in Figure 4: (a) moves rejected after evaluation at level 3 (solid line), (b) moves rejected after evaluation at level 2 (long-dash line), (c) moves rejected after evaluation at level 1 (short-dash line), and (d) moves evaluated with full potential (dotted line).
J. Phys. Chem., Vol. 100, No. 8, 1996 3269
Figure 7. Distribution of minimized energies as a function of rms deviation from the crystal structure, for 128 structures of 1CTF.
Figure 8. Distribution of minimized energies as a function of rms deviation from the crystal structure, for 128 structures of 1R69. Figure 6. Distribution of minimized energies as a function of rms deviation from the crystal structure, for 128 structures of 1MBO.
a factor of 5-10. The values of 〈Rg〉 are typically 10-15% greater than those of the crystal structures, taking into account the approximate packing of the calculated structures. This can be seen as a rough estimate of the effective balance between the attractive and repulsive parts of the potential. Some increase in size is an unavoidable consequence of the coarse-grained nature of the excluded volume, which makes it extremely difficult to achieve the precise packing of the native structure. The volume can be reduced by changing the cutoff distances, but at the expense of selectivity as the number of possible folds increases. There was some variation in this effect with the sequence, with 1R69 generating compact low-energy structures with a large overlap energy and few overlaps while 1CTF failed to minimize the energy even with a low overlap energy allowing more overlaps. The results are summarized in Table 1, where it can be seen that the average overlap energy per residue is roughly constant for each sequence. The energy coefficient of Rg itself was set to zero for all final results, with initial values in the range 5-10 used to accelerate the generation of compact structures and then annealed out. Typical distributions of the minimization results are shown in Figures 6-8 for ensembles of 128 structures. Results are shown for typical annealing runs with 3 × 105 trial moves per molecule. The temperature varied from an initial value on the
order of the standard distribution of the ensemble to a final value on the order of the difference between individual structures. Typical numbers were 20.0 and 1.0, respectively. These results show considerable variation in the performance of the model for the different sequences. The lowest root-mean-square (rms) deviation results are comparable to earlier results, but there is still little, if any, improvement in the correlation between the energy and the rms deviation which limits the predictive value of the results. In particular, the energy gap between the lowestenergy structure and the crystal structure varies considerably. In order to evaluate these results, it is useful to make a direct comparison to the previously published results in refs 2 and 3. Table 2 shows a comparison of the lowest-energy clusters of structures from the earlier work evaluated using the present potential function. These clusters represent groups of similar structures (as determined by their mutual rms deviations) corresponding to the distinct folded topologies produced in the simulations. In the case of 1R69, such a distinction could not be made objectively. Since the energy scale is different in the two cases, the values in Table 2 are scaled by the overall standard deviation of the distribution and shifted so that the average value is zero. The number of standard deviations separating the native energy from the calculated structures provides an index of the selectivity of the potential function. The results show that for 1MBO and 1CTF the energy gap is increased significantly, from a fraction of the standard distribution to a number on the order of two or three. In the case of
3270 J. Phys. Chem., Vol. 100, No. 8, 1996
Gunn
TABLE 2: Relative Energy Gaps of Previously Calculated Structures Screened with the Present Potential Functiona cluster
av rms deviation
1MBO(1) 1MBO(2) 1CTF(1) 1CTF(2) 1CTF(3) 1R69
7.74 12.41 4.78 9.33 8.73 8.85
av
old potentialb min
native
av
present potential min
native
-0.11 0.06 -0.30 -0.21 0.31 0.0
-2.98 -2.63 -3.74 -3.72 -3.58 -3.99
-3.01 -3.01 -4.25 -4.25 -4.25 0.80
-0.62 0.35 0.22 -0.60 0.36 0.0
-2.32 -1.87 -1.67 -2.25 -1.70 -3.31
-4.22 -4.22 -4.94 -4.94 -4.94 -2.42
a Energy values are presented in dimensionless units where the overall Average for each molecule is zero and the standard deviation of the distribution is one unit. b References 2 and 3.
1R69, the native energy was previously aboVe the average, and while there are still lower-energy structures with the new potential the native is more than two standard deviations below the averagesa qualitative improvement. In comparing the different topologies of 1MBO, the present potential shows an increased selectivity for the low-deviation structures, with a much larger gap between the average values although the distributions are still overlapping to a large extent. In the case of 1CTF, there is also a reordering of the clusters, although in this case disfavoring the low-deviation structures. This can be partially explained by the use of an explicit strand-pairing energy in ref 3 which artificially favored the native-like structures relative to the bare potential in that case. This suggests that some form of strand-pairing interaction may be required in addition to a generic contact potential in order to correctly distinguish the structures of β-sheet proteins. These results are similar to the preliminary results obtained for the first four-index tabular potential,3 which showed that such a potential could be very effective at screening the structures minimized with a cruder potential. The reestablishment of a substantial energy gap suggests that further minimization can continue to lead to improved structures. This was supported by a simulation carried out using an initial ensemble of 1MBO structures with very low rms deviations from the native (3.6 Å on average) which was generated by explicitly fitting the native structure to the reduced-model representation. Because these structures were determined without any potential function, they have many bad steric interactions which with earlier potentials caused the structures to be unstable, i.e., to unfold during minimization yielding a final distribution similar to that obtained from an extended initial configuration. With the present potential, however, many of these structures were stable, and the results of a minimization of these structures (with the same overlap energy as in Figure 6) is shown in Figure 9. This represents a considerable reduction in energy from the original fitted structures which had an average energy of +40 units. The resulting range of energies is comparable to that obtained from a random initial ensemble, and although many structures did unfold there are some which have low deviations and remained among the lowest in energy. This demonstrates that the desired solutions exist, that is, that there are nativelike structures within the model with competitive energies; however, there is no significant energy gap between these structures and many other alternatives. Because the volume of conformational space occupied by 14 Å rms structures is many orders of magnitude larger than that occupied by 4 Å rms structures, it is unlikely that low-deviation structures will be observed unless there is a large “basin of attraction” with an energy gradient leading to the native structure. The present results indicate, however, that in each case the apparent energy gap is quickly reduced or eliminated by direct minimization with the screening potential. Quantitatively, the best of the previously generated structures only account for 52%, 55%, and 75% (for 1MBO, 1CTF, and 1R69, respectively) of
Figure 9. Distribution of minimized energies as a function of rms deviation for 128 structures initially fit to the crystal structure of 1MBO, following energy minimization.
the present range of energies, with the difference in the case of 1MBO exactly equal to the relative lowering of the native energy with the new potential. However, the direct minimization shown here required roughly an order of magnitude less computer time than in the earlier work. The fact that the improved minimization leads to no real lowering of the average deviation has implications for the potential function. The initial screening results may be due more to a changing emphasis on different contributions to the energy than to a real increase in accuracy. As a result, the total energy is sensitive to small reorganizations of the structure which may not lead toward the correct structure in a useful way. A similar trend has also been observed16 with a number of other available potential functions. Structures generated with molecular dynamics were ranked differently when screened by different potentials, and in almost all cases the native structure was not correctly identified. Identification of the correct structure with a contact potential of this type requires that the combination of native contacts have the lowest energy; however, with a large number of possible contacts and a large number of statistically determined potential parameters, the accumulation of error in the total contact energy makes this unlikely.17 Most of the lowenergy structures generated here were found to differ significantly from one another, further suggesting that the increased detail of the potential leads to a greater number of distinct structures with comparable energies. The effectiveness of the potential might therefore be improved with some additional smoothing and increased emphasis on long-range interactions (physical or statistical). Better advantage could also be made of the multibody nature of the potential function, to better recognize favorable local conformations in addition to pairwise contacts as has been done with lattice models.18 This also implies that it can be very misleading to evaluate a potential function by screening a fixed set of structures. The
Minimizing Reduced-Model Proteins
J. Phys. Chem., Vol. 100, No. 8, 1996 3271
TABLE 3: Comparison of the Selectivity for the Native Structure at Each Level of Approximationa level ensemble of structures
0
1
2
3
minimized structuresb random compact structures
-2.03 -6.60
-1.09 -6.23
0.44 -2.19
0.16 -3.22
a The value of the native potential is given in terms of standard deviation units away from the mean of the ensemble. b As shown in Figure 6.
empirical contact potential upon which the current tabular potential is based was developed to identify correct motifs by threading sequences through the coordinates of all equal-length fragments in the database.11 In that case, the energies associated with a given sequence were found to have roughly a Gaussian distribution with the native structure several standard deviations lower. Roughly the same distribution was observed in running the present program without using the potential table, i.e., using only the overlap energy and a term proportional to the radius of gyration to generate otherwise random structures of the same secondary structure with the same compactness (16.0 average Rg, 1.3 average overlaps). The mean and standard deviation of this distribution were +9.4 and 28.7, respectively, compared to -113.4 and 32.8 for the results shown in Figure 6. The relative energy gaps for the minimized and random compact structures are shown in Table 3. The results for the random compact structures are comparable to the published results of threading; however, the results of minimization are qualitatively different, generating a nonoverlapping distribution at a much lower energy. This suggests that threading a sequence through the structures of other proteins is not far from a completely random sampling and therefore is a very poor indication of how the potential function will perform when it is actually minimized. Table 3 also shows the performance of the approximate potential functions with respect to the same distributions. As one would expect, the approximate potentials have a greatly reduced resolution, with the higher levels being completely unable to distinguish the native from the minimized structures. (This is expected, since their role is to screen unproductive trial moves, not the final distribution.) With respect to the other compact structures however, all levels of approximation show the same qualitative difference, with even the crudest approximation clearly being able to distinguish the native and minimized structures from the rest. The approximate potentials are thus capable of identifying completely incorrect folds even with a similar number of total contacts and compactness. The performance of the approximate potentials in distinguishing the two distributions further underlines the fact that the set of lowenergy structures, while still too large to be of predictive value, is many orders of magnitude smaller than the number of possible compact structures having the same secondary structure. V. Computational Performance The performance of the algorithm for different choices of the filtering criteria was investigated using a shortened test run of the 1MBO simulation. The parameters were chosen to correspond to an intermediate point in the simulation in order to make it as representative as possible. The energy scale, therefore, is not directly comparable to that of the previous section. Results are shown in Figure 4 for several choices of the adjustable cutoffs, with A being the most generous set and G the most strict. Two sets of results are shown, corresponding to different numbers of iterations. The first set consists of simulations all having the same total number of trial moves as the reference case where no screening was carried out. In all
TABLE 4: Comparison of Results for Several Test Runs of 128-Molecule Ensembles of 1MBO Using Different Cutoff Tolerances; Two Sets of Results Are Shown, Corresponding to Constant Iteration Count and Approximately Constant CPU Time cutoff four-index accepted av energyc level look-upsa movesb (arb units)
total CPU time
time per time per trial move accepted (ms) move (s)
10240
240
-32.6
30 782
23.48
1.00
A1 B1 C1 D1 E1 F1 G1
2636 936 641 498 383 232 103
226 148 141 104 82 68 35
-33.5 -23.6 -21.2 -17.6 -10.0 -5.6 5.3
14 309 6 363 5 551 4 590 3 650 3 416 2 580
10.92 4.85 4.24 3.50 2.78 2.61 2.17
0.50 0.34 0.31 0.35 0.35 0.39 0.58
A2 B2 C2 D2 E2 F2 G2
3441 2815 1922 1880 1781 924 357
303 435 399 375 344 228 107
-37.8 -41.1 -39.4 -36.1 -31.3 -27.3 -9.0
19 213 20 322 19 431 20 061 20 332 19 632 22 022
9.77 3.45 2.96 2.55 2.07 1.87 1.53
0.50 0.37 0.38 0.42 0.46 0.67 1.60
none
a This refers to the number of structures per molecule evaluated with the full potential. b Number per molecule. c See text for discussion of the energy scale used here.
cases there is some loss of accepted moves due to error in the screening, accompanied by a substantial reduction in CPU time. The efficacy of the minimization is judged by the final average energy obtained. The initial configuration in all cases was taken from the results of a short annealing calculation ending with the parameter values used here (overlap energy ) 3.0 and hydrophobic energy ) 5.0) and which had an average energy of +36.5 in dimensionless units. The results depend generally on the number of accepted moves and range from roughly half the reference reduction in energy to virtually the same result, obtained with choice A in less than half the CPU time. In practice, however, the simulations are run as long as reasonably possible, and the question is whether or not the potential can be minimized at all. A more appropriate test, therefore, is whether or not lower energies can be reached in a comparable time due to increased efficiency. In this case the results are clearly superior in all but the most extreme cases, with lower energies found in two-thirds the total time. The best results are obtained for intermediate values of the cutoffs, showing that there is a balance between the extremes of case A, where only a small increase in the number of iterations is possible, and case G, where the cutoffs are so restrictive as to impede the minimization. Except in the last two cases where the simulation appears to be getting trapped, there is relatively little change in the time per accepted move as the iteration counts are increased even though the overall acceptance rate is lower. The effect of choosing different levels of approximation was studied by repeating the same experiment with different segment lengths in the approximate potentials. The results are shown in Table 5 for several choices with the same number of iterations as in the reference in Table 4. The choice of levels used up to this point (9,3,1) is represented by case B1 described above. The cutoffs in the other cases were adjusted (to the extent possible) to yield a comparable number of accepted moves. The choice (9,3,1) requires the least amount of time to obtain this result, with fewer or more closely spaced approximations less efficient. The use of an even cruder level of approximation yields the most accepted moves, but at the expense of both speed and energy minimization. The use of only the single-index residue-residue potential as a filter results in a lower energy, but with a substantial increase in the total time.
3272 J. Phys. Chem., Vol. 100, No. 8, 1996
Gunn
TABLE 5: Comparison of Results for Several Choices of Levels of Approximate Potentials Used for Screening; Results Are for a 128-Molecule Ensemble of 1MBO with a Constant Number of Trial Moves coarse-graining levels 12 9 4
3 3 2 4 3
1 1 1 1 1 3 1
four-index look-upsa
accepted movesb
av energyc (arb units)
time per trial move (ms)
time per accepted move (s)
802 936 842 641 825 1346 706
213 148 200 152 162 151 147
-19.8 -23.6 -21.6 -13.1 -23.9 -22.6 -25.9
6.71 4.85 6.77 8.18 6.31 7.15 11.0
0.32 0.34 0.35 0.55 0.40 0.48 0.77
a This refers to the number of structures per molecule evaluated with the full potential. b Number per molecule. c See text for discussion of the energy scale used here.
TABLE 6: Comparison of CPU Times with and without Screening for Different Sequence Lengths molecule no. of residues CPU time (ms)/trial move
1MBO 1CTF 1R69 146 66 60
no screening 23.5 with screening 3.45 ratio 6.8 CPU time (s)/accepted move no screening 1.00 with screening 0.365 ratio 2.7
5.04 1.78 2.8 0.365 0.291 1.3
4.36 1.19 3.7 0.439 0.332 1.3
The scaling of the performance with sequence length was evaluated by carrying out a similar (although not as extensive) set of experiments for 1CTF and 1R69. In both cases, the choice of the (9,3,1) potential approximations was still found to be superior, suggesting that the efficiency depends on the spacing of the approximations, independent of the overall size. The fact that the approximate potentials were evaluated with fewer segments in the smaller molecules was reflected in the need for more generous cutoff criteria to obtain good results. The results of this comparison are summarized in Table 6 where the times per move are shown with and without screening. In all cases, lower energies were obtained in less time with the screening; however, the relative advantage was greater for 1MBO than for the smaller molecules. While the time per move in the unscreened case is very nearly proportional to N2, it is reduced to close to linear with the screening. Clearly one must be cautious drawing conclusions from three data points; however, these preliminary results suggest that greater improvements will be seen with larger molecules. These tests have not yet been undertaken since questions regarding the potential function make it difficult to interpret the results; however, such tests now appear to be feasible. VI. Conclusion The hierarchical table-look-up potential has been demonstrated for three examples of medium-sized proteins. The multiindex tabular potential table has a very general form which allows the correlations between the degrees of freedom to be explicitly taken into account. It has been shown that the table can be systematically generated by averaging over an all-atom potential function. The use of a series of further averaged tables as screening functions improves the overall efficiency of the program and makes the use of this form of potential practical for minimization. The increased specificity of the table thus leads to the generation of reasonable compact structures with reduced computational effort. These results show that the approximate versions of the potential used as screening functions, however, do not need to be particularly accurate in order to lead to a net advantage. All that is required is that the tests be relatively cheap and that the structures which are passed be “enriched” from the point of view of the full potential. The cutoff criteria can be continuously adjusted to optimize the overall performance. The method
described here is general and does not depend on the particular form of the data. While the effect is not as pronounced as in earlier work where an approximate functional form was known for the potential, it allows the same principle to be extended to more selective potential functions. The generality could be relaxed and thus the accuracy further enhanced by taking into account the particular structure of segments which are known to be rigid throughout the simulation. The results of the minimization show that the potential function used here is still not sufficient to reliably identify native-like structures for the molecules studied to date. The use of more efficient minimization techniques, however, provides the opportunity to further develop the potential function and apply more stringent tests of its selectivity for the native structure. The ability of a potential function to rank previously obtained structures does not appear to be a reliable indicator of how it performs when minimized directly. It is hoped that the current methods will encourage the development of other manybody potentials which can represent more of the details of amino acid structure while still taking advantage of the reduced-model representation. Acknowledgment. This work has been supported financially by the University of Montre´al, the Centre de Recherche en Calcul Applique´, FCAR (Que´bec), and NSERC (Canada), whose contributions are gratefully acknowledged. References and Notes (1) We ignore here the question of whether or not the native conformation really is the global minimum of a thermodynamic function and assume that an appropriate “effective” energy function can in principle be defined for which the observed structure is the global minimum. (2) Gunn, J. R.; Monge, A.; Friesner, R. A.; Marshall, C. H. J. Phys. Chem. 1994, 98, 702. (3) Monge, A.; Lathrop, E. A.; Gunn, J. R.; Shenkin, P. S.; Friesner, R. A. J. Mol. Biol. 1995, 247, 995. (4) Gunn, J. R.; Friesner, R. A. J. Comput. Chem., in press. (5) See for example: Multigrid Methods; Hackbusch, W., Trottenberg, U., Eds.; Springer-Verlag: Berlin, 1982. (6) Levitt, M.; Warshel, A. J. Mol. Biol. 1976, 106, 421-437. (7) Bashford, D.; Cohen, F. E.; Karplus, M.; Kuntz, I. D.; Weaver, D. L. Proteins: Struct., Funct., Genet. 1988, 4, 211-227. (8) It should be emphasized that this is a computational model and does not necessarily reflect actual folding mechanisms. (9) Casari, G.; Sippl, M. J. J. Mol. Biol. 1992, 224, 725. (10) Rooman, M. J.; Kocher, J. A.; Wodak, S. J. J. Mol. Biol. 1991, 221, 961. (11) Bryant, S. H.; Lawrence, C. E. Proteins: Struct., Funct., Genet. 1993, 16, 92. (12) Johnson, L. A.; Monge, A.; Friesner, R. A. J. Chem. Phys. 1992, 97, 9355. (13) Phillips, S. E. V. J. Mol. Biol. 1980, 142, 531. (14) Leijonmarck, M.; Liljas, A. J. Mol. Biol. 1987, 195, 555. (15) Mondragon, A.; Subbiah, S.; Alamo, S. C.; Drottar, M.; Harrison, S. C. J. Mol. Biol. 1989, 205, 189. (16) Wang, Y.; Zhang, H.; Li, W.; Scott, R. A. Proc. Natl. Acad. Sci. U.S.A. 1995, 92, 709. (17) Bryngelson, J. D. J. Chem. Phys. 1994, 100, 6038. (18) Kolinski, A.; Skolnick, J. Proteins 1994, 18, 338.
JP9524643
