HM Medical Clinic

 

Se359901368p

Global Optimization of Clusters, Crystals, and Biomolecules David J. Wales1* and Harold A. Scheraga2 crystals, and biomolecules. We focus on the Finding the optimal solution to a complex optimization problem is of great impor- PES, rather than the free energy, to avoid the tance in many fields, ranging from protein structure prediction to the design of additional complications arising from finite microprocessor circuitry. Some recent progress in finding the global minima of temperature. We also emphasize that, for potential energy functions is described, focusing on applications of the simple global optimization to succeed in predicting "basin-hopping" approach to atomic and molecular clusters and more complicated the properties of real systems, it is necessary hypersurface deformation techniques for crystals and biomolecules. These methods to have a sufficiently accurate potential ener- have produced promising results and should enable larger and more complex systems gy function. An efficient way to incorporate to be treated in the future.
more realistic potentials has been suggestedby Hartke (8).
The LJ model of inert gas clusters has he global optimization problem is a easily located because of the form of the been investigated intensively and provides a subject of intense current interest. Ap- potential energy landscape (2). For systems useful testing ground for putative global op- plications of obvious economic impor- of different sizes, locating the global mini- timization algorithms (14 ). In terms of the tance include traveling salesman–type prob- mum can be much harder and may require topology of the energy landscape and the lems and the design of microprocessor cir- sophisticated search routines. Similarly, a structure of stationary points, this potential cuitry. In the domain of atoms and molecules, random conformational search for the global also provides a useful model of noble gas discovering the lowest-energy isomer or crys- minimum of a typical protein would take an clusters (15 ). Furthermore, some of the non- tal structure for a system with a given com- unfeasibly long time. However, it is now icosahedral global minima first discovered position is frequently a goal. For example, it generally agreed that the search is not ran- for this potential have recently been identi- seems likely that the native structure of a dom, but instead the PES is biased toward the fied in nickel and gold clusters (16 ). A brief protein is often related to the global minimum global minimum, which results in efficient overview of global optimization strategies is of its potential energy surface (PES). Hence, relaxation (2, 3). Hence, the amino acid se- therefore first presented in the context of LJ considerable research efforts are being made quences of naturally occurring proteins are clusters (17 ). The prediction of crystal struc- to predict the three-dimensional structure of a presumed to have evolved to fold rapidly into tures is also a problem of current interest (18) protein solely from its amino acid sequence a unique native structure. This does not mean, and is discussed subsequently. Finally, we by computer simulation. Also, experimental however, that the corresponding computa- describe several strategies for tackling the microcalorimetry for free sodium clusters (1) tional search becomes trivial. Although there protein folding problem (6 ).
has revealed highly irregular thermodynamic is some evidence that it may be possible to properties as a function of size. To explain simulate protein folding by brute force mo- these results, the global potential energy min- lecular dynamics (4 ), much faster computers Most global minima for LJ clusters contain- imum, which is expected to be the favored will be needed before such tasks can be rou- ing fewer than 100 atoms are based on ico- structure for the low-temperature experiments, tinely undertaken.
sahedral packing (Fig. 1B). The exceptions, must first be identified. As these two exam- Numerous approaches to solving the glob- , serve as particularly inter- ples show, developing methods for treating a al optimization problem have been suggested esting test cases, because the corresponding diverse range of systems spanning the fields (5–13). One difficulty with the global optimi- energy landscapes consist of two families of of chemistry, biochemistry, and materials sci- zation literature is that it is spread over jour- structures (2, 19). At these sizes, the lowest- ence is important.
nals in many different fields; our citations energy minimum based on icosahedral pack- In treating any nontrivial global optimiza- reflect our knowledge of applications to prob- ing acts as a trap and is widely separated from tion problem, the principal difficulty arises lems in molecular science. Another problem the true global minimum (Fig. 1, A and C).
from the number of minima on the PES, is that detailed comparisons of different ap- Actually, it is quite easy to find the global which usually increases exponentially with proaches are rare, partly because gathering minima of these clusters by seeding the start- the size of the system. An example is the reliable statistics for nontrivial problems is ing geometry with a core of the appropriate cluster of 55 atoms interacting by a Lennard- time-consuming. In this article we focus on a morphology. Indeed, most of the lowest Jones (LJ) potential, LJ few particular strategies that we have found known minima up to LJ the number of minima (excluding permuta- to work best for problems involving clusters, tained by Northby by counting nearest-neigh- tional isomers) is at least 1010. Nevertheless,in this specific case the global minimum is Fig. 1. Global minima
of the LJ potential for (A) 38 atoms (truncated
1University Chemical Laboratories, Lensfield Road, octahedron), (B) 55 at-
Cambridge, CB2 1EW, UK. 2Baker Laboratory of oms (Mackay icosahe- Chemistry and Chemical Biology, Cornell University, dron), and (C) 75 atoms
Ithaca, NY 14853-1301, USA.
(Marks decahedron).
*To whom correspondence should be addressed. E- mail: [email protected] 27 AUGUST 1999 VOL 285 SCIENCE www.sciencemag.org bor interactions for icosahedral packing been used to locate global minima for LJ interval [0,1]. The temperature, T, becomes schemes (20). However, our focus here is on clusters of sizes up to 100 atoms, except for an adjustable parameter, but is not used for unbiased methods that may be transferable to those containing 75 to 77 atoms.
other systems.
A third approach uses genetic algo- In an application of the MCM procedure Simulated annealing (21) probably pro- rithms (31) that mimic the evolutionary to LJ clusters, all the lowest known minima vided the first generally applicable tech- process by evolving a "genetic code" based were located up to 110 atoms, and global nique for global optimization. In this ap- on the structure and using the concepts of minima for LJ , LJ proach the state of the system is followed fitness, mutation, and crossover. A variety were located for the first time in by simulation as the temperature is de- of different implementations have been unbiased searches (39). We have found that creased slowly from a high value, in the proposed (32). In discrete genetic algo- an unbiased genetic algorithm can also find hope that it will eventually come to rest at rithms, the conformational space is subject all the lowest LJ minima up to 110 atoms, the global potential energy minimum. A to a binary encoding corresponding to the if each member of the population is mini- new global minimum was located for LJ "genotype" (33). However, it is also possi- mized after every step (43). Use of the with this approach (22), but otherwise sim- ble to work directly in physical coordinate minimization step means that the same ulated annealing does not appear to have space, corresponding to the "phenotype" catchment basin landscape is searched.
been very successful for LJ clusters, even (34), and such a "continuous" genetic algo- This enabled Deaven and Ho (34 ) to find in more sophisticated forms (10, 23). The rithm located new global minima for LJ new global minima for LJ problem is that the free-energy global min- (35 ). The latter study also in- ilarly, the "two-level simulated annealing" imum can change at a temperature where cludes an implicit "catchment basin" trans- approach enabled Xue to find new global the energy barriers are too high for the system formation of the energy landscape, which and LJ ; this method to escape from a local minimum.
appears to be common to several of the An alternative approach is based upon more successful global optimization algo- without resetting the coordinates to those of "hypersurface deformation" where the func- rithms. This transformation can be separat- the current minimum (36 ). The "exponen- tional form of the potential energy is delib- ed from consideration of how the resulting tial tunneling" approach of Barro´n et al.
erately altered. Some transformations smooth surface is actually searched and is de- was also used to search the catchment basin the surface and reduce the number of min- scribed in the following section.
surface and to locate the truncated octahe- ima, thereby making the global optimiza- (37 ) (Fig. 1A). The same tion problem easier (10, 24 –28). However, authors also located some of the new global the global minimum of the smoothed sur- We now outline the particular hypersurface minima reported in (39) by a seeded search, face must then be mapped back to the real deformation that appears to be a constituent which again involves local minimization of surface, and this reversing procedure is the of several studies in which new global min- structures (38). Wolf and Landman suc- key problem associated with such ap- ima have been found for LJ clusters cessfully located the LJ global minima up proaches. It is now known that the global (35–39). This "basin-hopping" approach has using local minimization and a minimum can change rather dramatically proved to be useful for a range of atomic seeded genetic algorithm (44 ).
under some of the smoothing procedures and molecular clusters as well as biomol- To explain why the catchment basin (18, 29). Hence, it is necessary to couple ecules, and is easy to implement (40). The transformation is useful for global optimi- smoothing with an efficient local search functional form is explained in the Appen- zation, it is helpful to examine the thermo- procedure, which can be applied in map- dix, part A, and in Fig. 2. In this transfor- dynamics of the deformed landscape. For ping minima back from the deformed hy- mation, the potential energy for every point LJ clusters with nonicosahedral global min- persurface to the original one (18, 30). To in the catchment basin of each local mini- ima, Doye and Wales found that the catch- improve efficiency, more than one mini- mum becomes the energy of that minimum.
ment basin transformation broadens the oc- mum of the smoothed surface must be These catchment basins partition all of con- cupation probabilities of the different mor- tracked backwards (18, 30).
figuration space, and so the potential ener- phologies (19). On the original surface the gy can vary only in discrete steps when the overlap between these distributions is proach, the distance scaling method (DSM, geometry moves from one basin to another.
small, and so the probability of escape from see the section on biomolecules), has been The transformation must be combined with the large basin of icosahedral structures is applied with a reversing procedure that in- a search strategy, and Monte Carlo sam- rather low. Other groups have experiment- volves both structural perturbations and pling provides two possibilities if the struc- ed with different sampling distributions short molecular dynamics (MD) simula- ture is reset to that of the current local within the framework of simulated anneal- tions (27 ). A short molecular dynamics run minimum, or allowed to vary continuously.
ing, with the same intention of broadening was carried out for each minimum being We refer to search techniques coupled to transition regions (45 ).
tracked back at each step of the reversing the catchment basin transformation as ba- Some timings for the basic MCM pro- procedure, and another molecular dynamics sin-hopping. The "Monte Carlo plus energy cedure applied to LJ clusters are given in run was carried out for the resulting mini- minimization" (MCM) procedure of Li and Fig. 3 and Table 1 (40). There are two free mum on the undeformed surface. The de- Scheraga (41) corresponds to coordinate parameters: the temperature was fixed, and formation lowered barriers between mini- resetting and is generally found to be the the maximum step allowed in the perturba- ma, allowing the molecular dynamics pro- more effective approach (19, 39, 42). Steps tion of coordinates that precedes each min- cess to explore configurational space more are proposed by perturbing the latest set of imization was adjusted dynamically to give efficiently. Coupled in this way to molec- coordinates and carrying out a minimiza- a particular acceptance ratio. The mean ular dynamics, the DSM was used to locate tion from the resulting geometry. A step is number of basin-hopping steps and cpu the nonicosahedral global minimum for accepted if the energy of the new minimum, time required to find the global minimum at LJ ; however, it failed in some other cases , is lower than the starting point, E a fixed temperature, are shown in Fig. 3 for (27 ). Using local search and self-consistent is greater than E , then the step is LJ up to n ⫽ 74. The results are averages mapping of deformed and undeformed min- accepted if exp[(E )/kT] is greater over 100 different random starting points in ima, a deformation-based method has now than a random number drawn from the each case. In fact the optimal temperature www.sciencemag.org SCIENCE VOL 285 27 AUGUST 1999 increases somewhat with size, and so Fig. 3 effort would be greatly assisted if it were ture of benzene without prior assumption of does not correspond to the best possible possible to predict crystal structures from the space group.
parameterization for each cluster. Some the intermolecular potential alone. Unfor- Deformation procedures, initially adopt- more detailed statistics are presented in tunately, many compounds of interest are ed for treating biomolecules, have now Table 1 for selected sizes. These results polymorphic, exhibiting alternative pack- been applied to predict crystal structures illustrate how difficult it is to find the ing modes. This phenomenon causes par- without making use of ancillary informa- truncated octahedron for LJ , and how ticular problems for the pharmaceutical in- tion such as the space group (18). In this easy it is to find the Mackay icosahedron dustry, where polymorphic "impurities" work the diffusion equation method (DEM) can lead to undesirable physical properties, and DSM (Appendix, part B were used to It will certainly be possible to improve which, for example, led to litigation surround- predict the crystal structures of hexasulfur upon these results by exploring the trans- ing the drug Zantac (48). Polymorphism and benzene, which were treated as rigid formed landscape more efficiently, or per- may also provide a stringent test of empir- molecules (18). After fixing only the mo- haps by a quite different approach. One pos- ical intermolecular potentials.
lecular geometry and the interaction poten- sibility is to find pathways by means of true The application of global optimization tial, the unit cell dimensions, space groups, transition states (46 ). Another is to use mo- techniques to crystal structure prediction is and the number of molecules in the unit cell lecular dynamics to simulate the system's at an early stage of development. Some were all computed, and the experimental evolution between minimizations. The opti- studies have attempted to generate plausi- crystal structures were successfully locat- mal algorithm in any given case is almost ble starting points for energy minimization ed. For benzene, the calculation succeeded certainly system-dependent and may involve using common coordination geometries, even when the number of molecules in the a combination of different methods. Unfortu- most probable crystal symmetries, close- unit cell was set to twice the experimental nately, testing the efficiency of different al- packing arguments, and statistical correla- value, which made the global optimization gorithms for nontrivial problems, such as tions (49). Monte Carlo-simulated anneal- problem considerably harder.
and LJ , is rather time-consuming.
ing has also been considered (50) and mo- Genetic algorithms have also been applied lecular dynamics techniques have enabled to analyze powder diffraction data, which solvent and kinetic effects to be simulated does not require an intermolecular potential Crystal engineering is an important branch (51). Williams (52) has approached this at all. For example, Kariuki et al. correctly of materials science whose aim is to design global optimization problem by Monte predicted the previously unknown crystal solids with particular properties (47 ). This Carlo sampling to predict the crystal struc- structure of ortho-thymotic acid (53). Thisapproach could prove very useful when sin-gle crystals of sufficient size or quality for Fig. 2. Illustration of
routine structure solution are impossible to the E˜(X) energy land-
scape transformation A
in two dimensions. (A)
Original surface. (B)
Predicting the native structure of a protein Each local minimum from its amino acid sequence alone is an of E(X) corresponds to
area of intense current research. The poten- a plateau or catchment tial savings of experimental time and effort basin for E˜(X). The sur-
alone have stimulated a number of ap- faces are colored con- proaches: sequence-homology employs the sistently according to the energy. (D) View
known structures of sequences that are sim- of the transformed sur- ilar to the one in question (54 ), whereas face from above. (C)
Cut through the com- bined E˜(X) and E(X)
surfaces for the red boxed region shown in all the other panels.
The catchment basins can have complicated boundaries because of C
the finite step size used in the minimizations.
Fig. 3. Mean time required to find the global
minimum with the MCM procedure for LJn up to n ⫽ 74 (abscissa); the averages are over 100 random starting points in each case. The cpu time (seconds) in the lower (red) curve corresponds to a 250 MHz Sun Ultra II pro- cessor, whereas the number of basin-hopping steps in the (black) upper plot is dimension- less (ordinate).
27 AUGUST 1999 VOL 285 SCIENCE www.sciencemag.org threading uses energy (or energy-like) tion (CASP3) (63) proteins. Figure 4 illus- strated transferability between atomic and functions to compare the sequence with trates the quality of one of the blind pre- molecular clusters and has also yielded useful structural motifs from a database of known results for biomolecules (11). However, more structures (55 ). Some previous methods complicated methods such as conformational have combined the global optimization of a space annealing have been shown to perform potential energy function with constraints We have provided an overview of some se- better for SPA, ACB, the periplasmic protein based on secondary-structure prediction lected recent developments in the field of HDEA, and a bipartite transcriptional activa- global optimization as applied to clusters, tor MarA (59).
these have also been quite successful (56 ).
crystals, and biomolecules. For atomic and For biomolecules, development of better Here we focus on three global optimization molecular clusters the basin-hopping ap- potential functions is clearly a priority. Im- procedures, DEM (25), DSM (27, 57 ) and proach coupled to search strategies based on proved potential functions must better de- conformational space annealing (CSA) (58) Monte Carlo sampling or genetic algorithms scribe structures that exist in nature, where (see Appendix, part C ), which have recent- seems to work well. Unbiased algorithms can the global minimum is expected to be sepa- ly produced encouraging results without often treat systems with at least 100 atoms or rated by an energy gap from higher energy the use of any sequence or structure anal- molecules reliably, and we expect biased or structures. The quality of the potential is also ogies and databases (59).
seeded approaches to be useful for signifi- a serious issue in structure prediction for We have applied potential function defor- cantly larger systems.
mation to this global optimization problem, The basin-hopping approach has demon- With further improvements in both algo- with careful mapping between deformed andundeformed minima using local search tech-niques. The underlying principle is to locate Table 1. Mean time and number of basin-hopping steps taken to find the global minimum for selected
LJ clusters with the MCM procedure. The standard deviation is similar to the mean time in each case. The large regions of conformational space con- statistics were obtained for 1000 random starting points for each size with a fixed acceptance ratio of 0.5 taining low-energy minima by coupling them and temperature T* (reduced units). T* is the size-dependent optimal temperature that gives the shortest to some of the greatly reduced number of mean time for the given acceptance ratio; it was determined by varying T in steps of 0.1 and gathering minima on the highly deformed surface. The statistics for samples of 100 random starting points. The cpu times are for a 250 MHz Sun Ultra II DSM and the DEM have been implemented processor. rel., relative.
to carry out the deformation, each being ap-plied to a different part of the potential func- tion (see Appendix, part B).
The DSM has been applied to united- residue polyalanine chains with a length of up to 100 residues and to staphylococcal protein A (SPA) (57 ). It has successfully located low-energy structures of polyalanine chains, predicting that the most stable structure in the absence of solvent is a straight ␣-helix for up to 70 residues. For 70 – 80 residues the most stable form is bent in the middle of the ␣-he- lix and, from 80 residues upward, the most stable structure is a three-helix bundle. ForSPA, a minimum very close to the experi-mental structure has been located. These re-sults show that hypersurface deformation canbe applied to the conformational analysis ofglobular proteins.
The greatest advantage of the CSA method (see Appendix, part C ) is that itfinds distinct families of low-energy con-formations. It was successfully applied toall-atom polypeptide chains for the pen-tapeptide enkephalin (58) and to the 20-residue membrane-bound portion of melit-tin (60). An application using a united-residue representation (61) was also suc-cessful; the method located very lowenergy structures for the fragment of SPAconsisting of residues 10 through 55 andfor apo calbindin (ACB), with a root-mean-square deviation of 2.1 Å and 3.9 Å, re-spectively, from the ␣-carbon trace of theexperimental structure (62). CSA is also acore part of a recently developed hierarchi- Fig. 4. Superposition of the crystal (red) and predicted (yellow) structures of the CASP3 periplasmic
cal approach to protein-structure prediction protein HDEA. The C␣ atoms of the fragment included between residues Asp25 (D25) and Ile85 (I85) (59) first fully used on Critical Assessment were superposed. The root-mean-square deviation is 4.2 Å. Helices 3, 4, and 5 are indicated as H-3, of Techniques for Protein Structure Predic- H-4, and H-5, respectively (59).
www.sciencemag.org SCIENCE VOL 285 27 AUGUST 1999 rithms and potential energy functions, we 3. R. Zwanzig, A. Szabo, B. Bagchi, Proc. Natl. Acad. Sci. zation, and Machine Learning (Addison-Wesley, Read- expect to see solutions of previously intrac- U.S.A. 89, 20 (1992); P. E. Leopold, M. Montal, J. N.
ing, MA, 1989).
Onuchic, ibid. 89, 8721 (1992); M. H. Hao and H. A.
34. D. M. Deaven and K. M. Ho, Phys. Rev. Lett. 75, 288
table global optimization problems in many Scheraga, J. Phys. Chem. 98, 9882 (1994); A. Sali, E.
different fields (64 ).
Shakhnovich, M. Karplus, Nature 369, 248 (1994);
35. D. M. Deaven, N. Tit, J. R. Morris, and K. M. Ho, Chem. Note added in proof: Using a basin-hop- J. D. Bryngelson, J. N. Onuchic, N. D. Socci, P. G.
Phys. Lett. 256, 195 (1996).
Wolynes, Proteins 21, 167 (1995); J. P. K. Doye and
36. G. Xue, J. Glob. Opt. 1, 187 (1991).
ping approach, R. H. Leary has found a new D. J. Wales, J. Chem. Phys. 105, 8428 (1996).
37. C. Barro´n, S. Go´mez, D. Romero, App. Math. Lett. 9,
nonicosahedral global minimum of T␣ sym- 4. Y. Duan and P. A. Kollman, Science 282, 740 (1998).
(14 ).
5. C. D. Maranas and C. A. Floudas, J. Chem. Phys. 100,
38. C. Barro´n, S. Go´mez, D. Romero, ibid. 10, 25 (1997).
1247 (1994); M. Oresic and D. Shalloway, ibid. 101,
39. D. J. Wales and J. P. K. Doye, J. Phys. Chem. A 101,
9844 (1994); P. Amara and J. E. Straub, J. Phys. Chem. 5111 (1997).
99, 14840 (1995); D. Cvijovic and J. Klinowski, Sci-
40. An example computer code in Fortran can be down- A. A transformed energy landscape. The following trans- ence 267, 664 (1995).
formation of the energy landscape does not change the 6. H. A. Scheraga, Biophys. Chem. 59, 329 (1996).
41. Z. Li and H. A. Scheraga, J. Mol. Struct. (Theochem) relative energies of any minima: E˜(X) ⫽ min{E(X)}, where
7. J. A. Niesse and H. R. Mayne, J. Chem. Phys. 105,
179, 333 (1988).
X represents the three-dimensional vector of nuclear
4700 (1996).
42. R. P. White and H. R. Mayne, Chem. Phys. Lett. 289,
coordinates and "min" signifies that an energy minimiza- 8. B. Hartke, Chem. Phys. Lett. 258, 144 (1996).
463 (1998); D. J. Wales and M. P. Hodges, ibid. 286,
tion is carried out starting from X. The transformed
9. F.-M. Dittes, Phys. Rev. Lett. 76, 4651 (1996); J.
65 (1998); J. P. K. Doye and D. J. Wales, New J. Chem. energy, E˜(X), at any point, X, becomes the energy of the
Barhen, V. Protopopescu, D. Reister, Science 276,
22, 733 (1998).
structure obtained by minimization. Each local minimum 1094 (1997); P. K. Venkatesh, M. H. Cohen, R. W. Carr, 43. A. Markham and D. J. Wales, unpublished work.
is, therefore, surrounded by a catchment basin of constant A. M. Dean, Phys. Rev. E 55, 6219 (1997); A. F.
44. M. D. Wolf and U. Landman, J. Phys. Chem. A 102,
energy consisting of all the neighboring geometries from Stanton, R. E. Bleil, S. Kais, J. Comput. Chem. 18, 594
6129 (1998).
which that particular minimum is obtained. The overall (1997); G. A. Huber and J. A. McCammon, Phys. Rev. 45. C. Tsallis and D. A. Stariolo, Physica A 233, 395
energy landscape becomes a set of plateaus, one for each E 55, 4822 (1997); J. A. Niesse and H. R. Mayne,
(1996); I. Andricioaei and J. E. Straub, J. Chem. Phys. catchment basin (Fig. 2), but the energies of the local J. Comp. Chem. 18, 1233 (1997).
107, 9117 (1997); U. H. E. Hansmann and Y. Okamoto,
minima are unaffected by the transformation. Aside from 10. S. Schelstraete and H. Verschelde, J. Phys. Chem. A Phys. Rev. E 56, 2228 (1997).
removing all the transition state regions from the surface, 101, 310 (1997).
46. J. P. K. Doye and D. J. Wales, Z. Phys. D 40, 194
the catchment basin transformation also accelerates the 11. P. Derreumaux, J. Chem. Phys. 106, 5260 (1997).
(1997); N. Mousseau and G. T. Barkema, Phys. Rev. E dynamics, because the system can pass between basins all 12. T. Huber and W. F. van Gunsteren, J. Phys. Chem. A 57, 2419 (1998).
along their boundaries. Atoms can even pass through each 102, 5937 (1998); J. Schneider, I. Morgenstern, J. M.
47. G. R. Desiraju, Science 278, 404 (1997).
other without encountering prohibitive energy barriers.
Singer, Phys. Rev. E 58, 5085 (1998); W. Wenzel and
48. W. H. DeCamp, in Crystal Growth of Organic Mate- B. The diffusion equation and distance scaling methods. K. Hamacher, Phys. Rev. Lett. 82, 3003 (1999).
rials, A. S. Myerson, D. A. Green, P. Meenan, Eds. (ACS The DEM achieves a smoothing of a potential function 13. A number of books are also available, for example, R.
Proceedings Series, American Chemical Society, f(X) by transforming it into the function F(X, t), which is
Horst, P. M. Pardalos, and N. V. Thoai, Introduction to Washington, DC, 1996).
the solution of the diffusion equation with f(X) as the
Global Optimization (Kluwer Academic, Dordrecht, 49. A. Gavezzotti, Acc. Chem. Res. 27, 309 (1994); A.
starting condition and t (time) the deformation parame- Netherlands, 1995), and there is also a journal for Gavezzotti, J. Am. Chem. Soc. 113, 4622 (1991); B. P.
ter. The DSM, which is applicable to pairwise interactions, specialists, Journal of Global Optimization.
van Eijck, W. T. M. Mooij, J. Kroon, Acta Crystallogr. transforms the distance between the centers of interac- 14. D. J. Wales, J. P. K. Doye, A. Dullweber, F. Y. Naumkin, B51, 99 (1995); A. Gavezzotti and G. Filippini, J. Am.
tions according to the formula r˜ ⫽ (r ⫹ r0t)/(1 ⫹ t), The Cambridge Cluster Database. Available at http:// Chem. Soc. 118, 7153 (1996); C. B. Aakeroy, M.
where r0 is the equilibrium distance for a pairwise inter- Nieuwenhuyzen, S. L. Price, ibid. 120, 8986 (1998).
action. The whole procedure consists of macro-iterations 15. D. J. Wales, J. Am. Chem. Soc. 112, 7908 (1990).
50. H. Karfunkel and R. J. Gdanitz, J. Comput. Chem. 13,
in which the parameter t controls the deformation chang- 16. E. K. Parks et al., J. Chem. Phys. 107, 1861 (1997);
1171 (1992).
es between two extreme values, tmax and tmin (t ⫽ 0 R. L. Whetten et al., Adv. Mater. 8, 428 (1996).
51. A. Gavezzotti, Faraday Discuss. 106, 63 (1997).
corresponds to the original energy surface).
17. L. T. Wille, Annual Reviews of Computational Physics 52. D. E. Williams, Acta Crystallogr. A52, 326 (1996).
C. Conformational space annealing. The CSA method VII, D. Stauffer, Ed. (World Scientific, Singapore, in 53. B. M. Kariuki, H. Serrano-Gonza´lez, R. L. Johnston, (58, 60) searches the whole conformational space in its early K. D. M. Harris, Chem. Phys. Lett. 280, 189 (1997).
stages and then narrows the search to smaller regions with 18. R. J. Wawak, J. Pillardy, A. Liwo, K. D. Gibson, H. A.
54. P. K. Warme, F. A. Momany, S. V. Rumball, R. W.
low energy as the distance cut-off, Dcut, which defines the Scheraga, J. Phys. Chem. A 102, 2904 (1998).
Tuttle, H. A. Scheraga, Biochemistry 13, 768 (1974);
similarity of two conformations, is reduced. As for genetic 19. J. P. K. Doye and D. J. Wales, Phys. Rev. Lett. 80, 1357
T. A. Jones and S. Thirup, EMBO J. 5, 819 (1986); D. A.
algorithms (33), CSA starts with a preassigned number of Clark, J. Shirazi, C. J. Rawlings, Protein Eng. 4, 751
randomly generated and subsequently energy-minimized 20. J. A. Northby, J. Chem. Phys. 87, 6166 (1987).
(1991); M. J. Rooman and S. J. Wodak, Biochemistry conformations. This pool of conformations is called the 21. S. Kirkpatrick, C. D. Gelatt Jr., M. P. Vecchi, Science 31, 10239 (1992); M. S. Johnson, J. P. Overington, T. L.
bank. At the beginning, the bank is a sparse representation of 220, 671 (1983).
Blundell, J. Mol. Biol. 231, 735 (1993).
the entire conformational space. A number of dissimilar 22. L. T. Wille, Chem. Phys. Lett. 133, 405 (1987).
55. D. Fisher, D. Rice, J. U. Bowie, D. Eisenberg, FASEB J. conformations are then selected from the bank, excluding 23. J. Ma, D. Hsu, J. E. Straub, J. Chem. Phys. 99, 4024
10, 126 (1996); R. Goldstein, Z. A. Luthey-Schulten,
those that have already been used; they are called seeds.
(1993); J. Ma and J. E. Straub, ibid. 101, 533 (1994);
P. G. Wolynes, Proc. Natl. Acad. Sci. USA 89, 9029
Each seed conformation is modified by changing from one C. Tsoo and C. L. Brooks III, ibid. 101, 6405 (1994).
(1992) K. K. Koretke, Z. A. Luthey-Schulten, P. G.
variable to one-third of the total number of variables per- 24. F. H. Stillinger and T. A. Weber, J. Stat. Phys. 52, 1429
Wolynes, Protein Sci. 5, 1043 (1996).
taining to a contiguous portion of the chain; the new vari- 56. J. Skolnick, A. Kolin´ski, C. L. Brooks, A. Godzik, A. Rey, ables are selected from one of the remaining bank confor- 25. L. Piela, J. Kostrowicki, H. A. Scheraga, J. Phys. Chem. Curr. Biol. 3, 414 (1993); A. Kolin´ski and J. Skolnick,
mations rather than being picked at random. Each confor- 93, 3339 (1989).
Proteins Struct. Funct. Genet. 18, 338 (1994); J.
mation is energy minimized to give a trial conformation.
26. J. Kostrowicki, L. Piela, B. J. Cherayil, H. A. Scheraga, Skolnick, A. Kolin´ski, A. R. Ortiz, J. Mol. Biol. 265, 217
For each trial conformation, ␣, the closest conformation A ibid. 95, 4113 (1991); F. H. Stillinger and D. K. Still-
(1997); B. A. Reva, A. V. Finkelstein, J. Skolnick, Fold. from the bank (in terms of distance D␣A) is determined. If inger, ibid. 93, 6106 (1990); T. Head-Gordon, F. H.
Des. 3, 141 (1998).
D␣A ⬍ Dcut (Dcut being the current cut-off criterion), ␣ is Stillinger, J. Arrecis, Proc. Natl. Acad. Sci. USA 88,
57. J. Pillardy, A. Liwo, M. Groth, H. A. Scheraga, J. Phys. considered similar to A; in this case ␣ replaces A in the bank 11076 (1991); R. J. Wawak, M. M. Wimmer, H. A.
Chem. 103, 7353 (1999).
if it is also lower in energy. If ␣ is not similar to A, but its Scheraga, J. Phys. Chem. 96, 5138 (1992); H. A.
58. J. Lee, H. A. Scheraga, S. Rackovsky, J. Comput. Chem. energy is lower than that of the highest-energy conforma- Scheraga, Int. J. Quant. Chem. 42, 1529 (1992); J.
18, 1222 (1997).
tion in the bank, B, ␣ replaces B. If neither of the above Pillardy, K. A. Olszewski, L. Piela, J. Mol. Struct. 59. A. Liwo, J. Lee, D. R. Ripoll, J. Pillardy, H. A. Scheraga, conditions holds, ␣ is rejected. The narrowing of the search (Theochem) 270, 277 (1992).
Proc. Natl. Acad. Sci. U.S.A. 96, 5482 (1999).
regions is accomplished by setting Dcut to a large value 27. J. Pillardy and L. Piela, J. Phys. Chem. 99, 11805
60. J. Lee et al., Biopolymers 46, 103 (1998).
initially (usually one-half of the average pair distance in the 61. A. Liwo et al., J. Comput. Chem. 19, 259 (1998).
bank) and gradually diminishing it as the search progresses.
28. J. Pillardy and L. Piela, J. Comp. Chem. 18, 2040
62. J. Lee, A. Liwo, H. A. Scheraga, Proc. Natl. Acad. Sci. Special attention is paid to selecting seeds that are far from (1997); M. A. Moret, P. G. Pascutti, P. M. Bisch, K. C.
U.S.A. 96, 2025 (1999).
each other. One round of the procedure is completed when Mundim, ibid. 19, 647 (1998).
63. Third Community Wide Experiment on the Critical there is no seed to select, that is, all conformations from the 29. J. P. K. Doye, D. J. Wales, R. S. Berry, J. Chem. Phys. Assessment of Techniques for Protein Structure Pre- bank have already been used. The round is repeated a 103, 4234 (1995).
diction (CASP3). Available at http://predictioncenter.
predetermined number of times.
30. R. J. Wawak, K. D. Gibson, A. Liwo, H. A. Scheraga, Proc. Natl. Acad. Sci. USA 93, 1743 (1996).
64. D.J.W. is grateful to J. Doye for his comments on this References and Notes
31. J. H. Holland, Adaptation in Natural and Artificial manuscript and to the Royal Society and the Engi- 1. M. Schmidt, R. Kusche, W. Kronmu¨ller, B. von Issen- Systems (Univ. of Michigan Press, Ann Arbor, 1975).
neering and Physical Sciences Research Council for dorff, H. Haberland, Phys. Rev. Lett. 79, 99 (1997).
32. A. A. Rabow and H. A. Scheraga, Protein Sci. 5, 1800
financial support. H.A.S. is grateful to J. Pillardy for 2. D. J. Wales, M. A. Miller, T. R. Walsh, Nature 394, 758
help in writing parts of this manuscript and to the 33. D. E. Goldberg, Genetic Algorithms in Search, Optimi- NIH and the NSF for financial support.
27 AUGUST 1999 VOL 285 SCIENCE www.sciencemag.org

Source: http://www-wales.ch.cam.ac.uk/pdf/SCIENCE.285.1368.1999.pdf