J . Phys. Chem. 1992, 96,4612-4616
4612
Mean Field Theory as a Tool for Intramolecular Conformational Optimization. 1. Tests on Terminally-Blocked Alanine and Met-Enkephalin K.A. Olszewski,*vt L.Piela,**'and H.A. Scheraga*q* Quantum Chemistry Laboratory, Department of Chemistry, Warsaw University, Pasteura 1 , 02-093 Warsaw, Poland, and Baker Laboratory of Chemistry, Cornell University, Ithaca, New York 14853-1 301 (Received: September 16, 1991)
A self-consistent multitorsional field (SCMTF) method of global optimization of the conformational energy of an oligopeptide is proposed. The method is based on the idea that the maximum of the square of the ground-state wave function is very often close to the global minimum of the potential energy and may be relatively unaffected by the complexity of the potential, e.g., its many minima. Using a mean field approximation, a set of coupled Schrijdinger equations for the motion of the nuclei is solved iteratively in a dihedral-angle space, each equation depending on a single dihedral angle. To obtain the effective potential energy, and therefore the Hamiltonian, for the particular dihedral angle, a Monte Carlo averaging is carried out over all other dihedral angles. Thus, the wave function and therefore the probability distribution of a particular dihedral angle depends on the mean field of the rest of the molecule. The method was found to converge rapidly when applied to searches for the global minimum of the potential energy of terminally-blocked alanine. In the case of the pentapeptide Met-enkephalin, convergence was slightly slower, but it also led to the global minimum. Because the number of evaluations of the potential energy increases approximately linearly with the number of dihedral angles (not counting the cost to compute the potential energy), the method is expected to be applicable to larger molecules.
1. Introduction
close to the global-minimum position of the potential energy. Thus, the wave function can identify the global minimum even though The search for the global minimum of the conformational energy of a macromolecule encounters the multiple-minima problem, Le., the existence of numerous minima in the multivariable energy hyperspace.' While algorithms are available for (1) Gibson, K. D.; Scheraga, H. A. In Structure and Expression, Volume minimizing a function of many variables, none exist for passing 1, From Proteim to Ribosomes; Sarma, M. H., Sarma, R. H., Eds.; Adenine Press: Guilderland, NY, 1988; p 67. effectively from one local minimum, over a potential barrier, to (2) Simon, I.; NBmethy, G.; Scheraga, H. A. Macromolecules 1978, 11, the next local minimum and ultimately to the global minimum, 797. or even for identifying the global minimum when it has been (3) Visquez, M.; Scheraga, H. A. Biopolymers 1985, 24, 1437. attained. Various procedures have been developed to surmount (4) Visquez, M.; Scheraga, H. A. J. Biomol. S t r u t . Dyn. 1988,5, 705. thii problem. The following are the most promising ones that have ( 5 ) Visquez, M.;Scheraga, H. A. J. Biomol. S t r u t . Dyn. 1988, 5, 757. been developed in our laboratory, and applied successfully, in most (6) Piela, L.;Scheraga, H. A. Biopolymers 1987, 26, S33. cases to small oligopeptides containing up to 20 amino acid res(7) Li, Z.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1987,84,6611. (8) Li, Z.; Scheraga, H. A. J . Mol. Struct. (THEOCHEM) 1988, 179, idues: (1) b u i l d ~ p , ~(2) - ~ optimization of electrostatics (self333. consistent electrostatic field, SCEF),6 (3) Monte Carlo plus energy (9) Ripll, D. R.; Scheraga, H. A. Biopolymers 1988, 27, 1283. minimi~ation,~** ( 4 ) electrostatically driven Monte Carlo (EDMC),+I' ( 5 ) adaptive importance sampling Monte C a r l ~ , ~ ~ - ' ~ (10) Ripoll, D. R.; Scheraga, H. A. J . Protein Chem. 1989, 8, 263. (11) R i p l l , D. R.; Scheraga, H. A. Biopolymers 1990, 30, 165. ( 6 ) relaxation of dimen~ionality,'~.~~ (7)pattern recognition im(12) Paine, G. H.; Scheraga, H. A. Biopolymers 1985, 24, 1391. portance sampling minimization (PRISM),I7-l9 (8) diffusion (13) Paine, G. H.; Scheraga, H. A. Biopolymers 1986, 25, 1547. equation method,20g21and (9) use of h o m o l ~ g y . ~It~is* still ~ ~ an (14) Paine, G. H.; Scheraga, H. A. Biopolymers 1987, 26, 1125. open question, however, as to whether these methods can be scaled (15) Purisima, E. 0.;Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A.1986, 83, 2782. up to treat proteins containing 100-200 amino acid residues. (16) Purisima, E. 0.;Scheraga, H. A. J. Mol. Biol. 1987, 196, 697. In this paper, we propose a self-consistent multitorsional field (17) Lambert, M. H.; Scheraga, H. A. J . Compur. Chem. 1989,10,770. (SCMTF) method, based on solving the Schrcidinger equation for (18) Lambert, M. H.; Scheraga, H. A. J . Comput. Chem. 1989,10, 798. internal rotational motions (variations of dihedral angles) in the (19) Lambert, M. H.;Scheraga, H. A. J . Compur. Chem. 1989,10,817. molecule under study. Bond lengths and bond angles are main(20) Piela, L.; Kostrowicki, J.; Scheraga, H. A. J . Phys. Chem. 1989, 93, tained fixed as in the empirical conformational energy program 3339. for peptides (ECEPP) a l g ~ r i t h m , * and ~ - ~a~ mean field approx(21) Kostrowicki, J.; Piela, L.; Cherayil, B. J.; Scheraga, H. A. J . Phys. imation is used; consequently, one has to solve a set of coupled Chem. 1991, 95, 4113. (22) Warme, P. K.; Momany, F. A,; Rumball, S. V.; Tuttle, R. W.; one-dimensional equations, the variable in each equation being Scheraga, H. A. Biochemistry 1974, 13, 768. a single dihedral angle. The procedure is repeated iteratively until (23) Acharya, K. R.: Stuart, D. I.; Phillips, D. C.; Scheraga, H. A. J. self-consistency is achieved for all the dihedral angles of the Protein Chem. 1990, 9, 549. molecule. No solvent effects have been taken into account ex(24) Momany, F. A,; McGuire, R. F.; Burgess, A. W.; Scheraga, H. A. plicitly in the present paper except for some screening of the J. Phys. Chem. 1975, 79, 2361. electrostatic potential in the ECEPP algorithm (dielectric constant (25) Ntmethy, G.; Pottle, M. S.; Scheraga, H. A. J . Phys. Chem. 1983, 87, 1883. equal to 4). Solvent effects can be included relatively easily by (26) Sippl, M. J.: Ntmethy, G.; Scheraga, H. A. J . Phys. Chem. 1984,88, using a method for calculating the water-exposed surface in the 623 1. molecule.27 (27) Perrot, G.; Cheng, B.; Gibson, K. D.; Vila, J.; Palmer, K. A,; Nayeem, The general philosophy of our procedure is similar to that of A.; Maigret, B.; Scheraga, H. A. J. Comput. Chem. 1992, 13, I . Ebeling et al.,2s RujaqZ9and S o m ~ r j a iwho , ~ ~used ~ ~ the ~ well(28) Ebeling, W.; Engel, A,; Esser, B.; Feistel, R. J . Stat. Phys. 1984, 37, 369. known ruleg2that the position of the maximum of the square of (29) Rujin, P. Z . Phys. B: Condens. Matter 1988, 73, 391. the ground state wave function for any kind of system may be
* Cornell University.
'Warsaw University.
*Towhom requests for reprints should be addressed at Cornell University.
(30) Somorjai, R. L. J . Phys. Chem. 1991, 95, 4141. (31) Sylvain, M.; Somorjai, R. L. J. Phys. Chem. 1991, 95, 4147. (32) Landau, L. D.; Lifshitz, E. M. Quantum Mechanics; Pergamon Press: New York, 1958; Chapter VII.
0022-365419212096-4672%03.00/0 0 1992 American Chemical Society
The Journal of Physical Chemistry, Vol. 96, NO. 11, 1992 4673
Intramolecular Conformational Optimization the potential energy may have many local minima. In order to increase the volume of conformational space explored, the authors regulate the partitioning between the kinetic energy and the potential energy, a procedure which we do not always find to be necessary. Ebeling et a1.,2* as well as R ~ j a n , 2did ~ not relate their global optimization method to any molecular problem. Somorjai's method and our procedure are designed to deal with molecules. There are, however, differences between our procedure and that of Somorjai. First, we solve a multidimensional intramolecular rotational problem, whereas Somorjai thus far has applied his procedure to one-dimensional mathematical problems. Second, our method is a self-consistent mean field theory, whereas the potential energy in Somorjai's approach is fixed throughout his procedure. More importantly, if Somorjai's global minimization procedure were applied to a multidimensional problem, one would still be left with the necessity to determine the location of the global maximum of the total wave function, a problem of complexity similar to the original one. Also in Somorjai's procedure, applied thus far to one-dimensional problems, use is made of Gaussians to represent the ground-state wave function; this requires that the Gaussians (or any other chosen functions) form a complete set. For real problems, however, the basis set has to be finite, but as large as possible. For a molecule, the number of Gaussians required is proportional to L3M,where L is the maximum size of the molecule in angstroms (assuming that 1 A is a reasonable uncertainity for the position of any atom) and M is the number of atoms. Since the diagonalization of the Hamiltonian matrix required in his method, therefore, varies as L9M-', it can be seen that his method is not applicable at present to molecular problems. The use of the one-center expansion proposed by Somorjai does not alleviate this problem because a large distance between the expansion center and the atomic nuclei leads to an exponential increase in the required number of expansion functions.
2. Self-consistent Multitorsional Field Method Consider a molecule composed of M atoms with masses m,, (n = 1, ...,Mj. The Hamiltonian operator for the motion of the nuclei is
ri=i.+t
(1)
where the circumflex over the symbols means that the corresponding quantity is an operator. In the Cartesian coordinate representation, the kinetic energy operator T i s given by
where A,, is the Laplacian operator for the nth particle and v i s the potential energy operator. Since bond lengths and bond angles usually vary much less than dihedral angles, it is reasonable to approximate them as constant at their equilibrium values. Therefore, the only remaining variables are the dihedral angles 0 = (01, ..., In order to solve the Schriidinger equation in dihedral-angle space, Le., the space of the molecular conformations with frozen bond lengths and bond angles, we have to transform the Hamiltonian of eq 1 appropriately. In dihedral-angle space, the potential energy operator has to be expressed in terms of dihedral angles, and the kinetic energy operator takes the form33*34 &)e
(3)
ai
ai
where is the momentum operator conjugated with (the position operator) and g = det (C) is, in general, a function of 8 defined by the matrix C-l given by (4) (33) GO, N.; Scheraga, H. A. Macromolecules 1976, 9, 535. (34) Bunker, P. R. Molecular Symmetry and Spectroscopy; Academic Press: New York, 1979; p 121.
where r,, is the Cartesian coordinate vector representing the position of the nth atom in the molecule. The procedure of Noguti and G635is used to compute G-I, which is then inverted to produce G. The nondiagonal elements of C were found here to be at least 1 order of magnitude smaller than the diagonal ones. Therefore, in order to demuple the variables 8, in the kinetic energy operator, we make an approximation that the nondiagonal elements of G are equal to zero, and that G is independent of 8, Le. Gij N (1
(5)
where Ii is a moment of inertia and 6, is the Kronecker delta. With this approximation we obtain a model Hamiltonian
which is
(7) in &representation. In a first approximation, as in the Hartree method in quantum theory,36we may treat the dihedral angles as independent of each other (mean field approximation) as follows. We assume that the solution of the corresponding eigenproblem may be approximated by the Hartree-like product of the normalized one-angle functions 4,(ei)
weI, e2, ..., ON) = fi4,(ei) I i=
(8)
According to the variational principle E = (*l@&l\k) L Eo
(9) where Eo is the ground-state energy of the molecule. Therefore, by minimizing E by varying 4i, i = 1, ..., N, one obtains a set of N-coupled one-dimensional equations,36viz.
H& = &+
i = 1,
..., N
(10) where 4p means the kith excited state (ki = 0 is the ground state; those states are needed for insertion in eq 8 in order to obtain the total ground-state wave function \to)and
rii = Pi + pff(ei)
(11)
where the probability density distribution pp 14fI2. It can be seen from eq 13 that the effective potential Vy"(0,) depends on the mean field created by averaging over the other dihedral angles, OI (1 # i). In order to calculate Vy", the integrals appearing in eq 13 have been approximated with a Monte Carlo procedure as A =
1 M qff(e,)-E v(e;.',..., e?;, ei, e$"', ..., em,') MEm
(14)
where the summation extends over MEtrial points in the N - 1 dimensional space of all dihedral angles Or except Oi. First, each point 8" in the space is selected randomly according to some preassumed onedimensional distribution p f , 1 = 1,2, ,..,N. Then, whenever a Om is chosen, the potential energy V(O)is minimized with respect to all 8's. The minimization gives the new P*,which is then used in eq 14. Solution of eqs 10 gives a new set of 47 and, therefore, a new set of pp. The procedure is repeated iteratively until self-consistency of p: distributions is achieved. This approach is by far less expensive than the full scanning of the (35) Noguti, T.; G6,N. J . Phys. SOC.Jpn. 1983, 52, 3283. (36) Landau, L. D.; Lifshitz, E. M. Quantum Mechanics; Pergamon Ress: New York, 1958; pp 232-235.
4674 The Journal of Physical Chemistry, Vol. 96, No. 1 1 , 1992
Olszewski et al.
has the meaning of the mean value of the potential energy. In summary, the algorithm consists of the following steps: (1) Assume starting probability distributions for each Oi. (2) Determine li of eq 5 from the procedure described below, which involves generating conformations from th,e robability distributions pp(Oi), i = 1 , ...,N . (3) Construct &Bi) by using eq 14 for each dihedral angle Bi. The 0" are generated from the probability distributions pp(e,), and Om* are obtained by minimization. Keep a record of the lowest potential energy conformation (lowest V) that appeared in the minimization, if it is lower than anyone found until now. (4) Solve the set of s 10 and then compute pp = l$Pl2. (5) Go to step 2, unless the 'pi 8 distributions no longer change. (6) The resulting geometry is that of the lowest energy calculated in step 3 during the whole procedure. Concerning step 1 the starting density distribution for the first iteration was chosen as the Dirac 6 function ,,jni'(e)
Dihedral A n g l e (degrees)
Figure 1. Effective ECEPP potentials fo; Met-enkephalin (energy in arbitrary units). (A) A typical plot of Vff(ei)(for Oi = 4 of Gly-3) computed for Mc = 10 (a) and Mc = 100 (b); see eq 15. Both plots are practically the sate and superimpose very well. (B) The largest discomputed with M c= 10 (a) and Mc = 100 (b), crepancy between which occurred for Bi = J/ of Met-5. Even in this case, the global-minimum well for this dihedral angle (around 160') is very accurately represented. Similar plots for terminally-blocked alanine show that the results are similarly nearly independent of M'. These results for N = 3 and N = 10 suggest that Mc is, to a good approximation, independent of N .
vff(Oi)
conformational space, because fiff(Bi) turned out to be quite insensitive to variations of Mc (see Figure 1). The idea behind eq 14 is the following. The effective potential fiff(fli)is created by averaging over the positions of the other dihedral angles with their probability distribution pp, 1 # i. The maxima of py correspond to the minima of the potential energy V(8),that is attained in eq 14 by using Om*. Thus, by staying away from very high values of V(O),we get rid of possible divergences that appear in mean field theories when evaluating the integral of eq 13. This problem arises when one uses a potential containing Lennard-Jones repulsive terms or any other hard-core rep~lsion.~' A somewhat similar method for eliminating hard-core regions from the potential is discussed in refs 38 and 39. The reason for the divergence is nonphysical, because the wave functions produced by the mean field equations exaggerate the importance of the hard-core repulsion. When the repulsive potential increases faster than 1/r around r = 0, it produces infinities in the integrals of eq 13, despite the fact that the probability of encountering such a situation is very small because the corresponding energy is extremely high. Once the effective potential energy is obtained, the mean value of the energy is calculated as
where $ = cy - (dpIeffI&)is interpreted as the mean value of the kinetic energy related to the ith dihedral angle, and (37) Thouless, D. J. The Quantum Mechanics of Many-Body Systems; Academic Press: New York, 1961; Chapter IV. (38) Moszkowski, S . A.; Scott, 9. L. Ann. Phys. 1960, 11, 65. (39) Moszkowski,S. A. Phys. Reu. 1965, B140, 283.
=
s(e - @nit)
(17)
where Bjni' ( i = 1, ..., N) is an arbitrary initial dihedral angle. Finally, after convergence was attained, the results did not depend on the starting distributions pinit. As to step 2, in order to determine li of eq 5 , the following procedure has been adopted. The matrix G-'was calculated from eq 4 and averaged over M' trial conformations (typically from 10 to 100) in order to eliminate the dependence on 8. Then the matrix was inverted, and the diagonal elements of G were taken as 1 /Ii; see eq 5 . Concerning step 3, a number (M') of trial conformations were generated and their energies were minimized locally; only a few of them (Me),those correspon9ing to the lowest energy obtained, The number of evaluations of were chosen to approximate the potential energy is K P N , where Kg = 640 was assumed in the present work. The last part of the computations consumes almost 100%of the total computer time; therefore, assuming that Mc is independent of N (see Figure l), the computational effort of the method increases linearly with the number of the dihedral angles involved in the calculations. The trial conformations were generated with the von Neumann algorithm.40 The ECEPP/2 p o t e ~ ~ t i a was I ~ ~assumed -~~ for V(0) and, therefore, the total computational effort scales as N4because the calculation of the ECEPP potential itself scales with N as M . As to step 4, the differential equations (10) were solved by using the renormalized Numerov-Cooley p r o c e d ~ r ewith ~ ' ~ the ~ ~ step size equal to 2s/K,. Then, the probability densities py were constructed from the solutions of eq 10 and used to generate another set of the M' trial conformations and so on. As expected in our calculations, the ground-state yave function $7 indicated the global minimum potential well of unless it was very narrow. However, it turned out that in some cases other low-lying minima were of significant importance, when searching for the global minimum of V. Therefore, to enhance the exploration of the conformational space, our procedure has been modified by applying the Boltzmann-averaged probability density pi( Rei) instead of pp(8,)
v'".
$7
vff((ei)
pi( T;B,) =
Z-1
k,=O
e-8(r,*r.0i)ld:i(ei)(2
(18)
where /3 = l/kTand Z'is a normalization factor. This procedure, involving a temperature-dependent density and the excited states, is often used in the literature, e.g., ref 43. Thus, the temperature T becomes a parameter of the procedure. At T = 0, one has pi(RBi)= pp(Bi), as before. Usually, a few consecutive runs, hereafter called macroiterations, are carried out, each with a (40) von Neumann, J . Natl. Bur. Stand. Appl. Math. 1951, Ser. 12, 36. (41) Cooley, J. W . Math. Comput. 1961, 15, 363. (42) Johnson, B. R. J . Chem. Phys. 1977, 67, 4086. (43) Kirzhnits, D.A. Field Theoretical Method in Many-Body Systems; Pergamon Press: New York, 1967; Chapter 5.
The Journal of Physical Chemistry, Vol. 96, NO. 11, 1992 4615
Intramolecular Conformational Optimization different value (but fixed during the macroiteration) of the parameter. Three different criteria for numerical convergence of the iterative procedure were tested. The criteria are related to the following sets of dihedral angles, respectively: ejnin = {e:qff(e)= min) (19) i.e., the 0 corresponding to the minimum value of the mean field potential. = (O:pp(t9) = max)
(20)
Le., the t9 corresponding to the maximum value of the probability density.
0 sin e7
4
= arctan
CL :I
COS
e7
i.e., the mean value of 0 extracted from the probability distribution pi( T;B) (after averaging the direction of the unit vector). The procedure was terminated when the values of the dihedral angles in any chosen set were close in succeeding iterations. This was an indication of the convergence of the corresponding density distributions pi( TOi). Finally, we comment on the question as to whether the maximum value of 1O0l2indicates the position corresponding to the global minimum of the potential energy or the position preferred by the free energy. If the potential well of the global minimum is not too narrow, then the maximum value of 1q0l2 indicates the global minimum of the potential energy function. However, if the potential well of the global minimum is narrow, then the ground-state wave function localizes instead in the broader potential well of higher energy, even though some excited-state wave functions may be partially concentrated in the vicinity of the global-minimum potential well. In the temperature-dependent Hartree approximation, eq 18, one directly takes into account not only the ground but also the excited states (unless T = 0) by constructing the Boltzmann-averaged density pi This means that pi will have maxima where the conformations preferred by the free energy are located. Since the conformations of eq 14 are generated according to the density p i , then the Vff(Oi)will qualitatively have the behavior of a free energy rather than an energy. Therefore, the ground-state solution of the Schriidinger equation (10) will be localized within a low (or the lowest) energy well, as indicated by free energy. Similar questions are discussed in ref 30.
3. Numerical Results and Discussion Two molecules were chosen in order to test the method: Nacetyl-N'-methylalanineamide (terminally-blocked alanine) (NANMA) and Met-enkephalin (ENK). The latter molecule is the opioid pentapeptide H-Tyr-Gly-Gly-Phe-Met-OH. Conformations related to the global-minimum ECEPP/2 energies of NANMAM and ENK7J6are known. For NANMA, the dihedral angles 4, $, and x were varied; the other dihedral angles (the two w's and the dihedral angles of the terminal CH3groups) were fmed at the global minimum values.M Calculations for enkephalin were carried out with 10 variable dihedral angles, i.e., 4i and $i for (i = 1, ..., 5 ) . The other dihedral angles, viz., wi,i = 1, ..., 5 , and the side-chain dihedral angles were frozen at the values corresponding to the global minimuma7J6 In both examples, p of eq 18 was used, and T was taken as 293 and 0 K. Terminally-BlockedAlanine. In the case of terminally-blocked alanine, the SCMTF method was convergent from any point of the conformational space. The 4,$-map of NANMA has eight local minima, and starts from any of these minima led to the global minimum in the first or second iteration. Table I contains an example of results of SCMTF iterations carried out a t T = 293 K (see eq 18) with Mc equal to 10. The case of Mc = 100 was also tested and did not appear to introduce any changes. Table (44) Viquez, M.;Nlmethy, G.; Scheraga, H. A. Macromolecules 1983,
16, 1043.
TABLE I: Application Alanine* iteration no. 1 2 3 4 5 6
of the SCMTF Method to Terminally-Blocked -
tr
0.21 0.42 0.42 0.42 0.42 0.42
t&
t"
0.19 0.60 0.60 0.60 0.60 0.60
0.21 2.41 2.41 2.40 2.41 2.41
P -4.30 -3.37 -3.36 -3.37 -3.36 -3.36
a T was taken as 293 K with Me = 10. The values of denote the mean kinetic energy corresponding to the ith dihedral angle, and is the mean value of the potential energy (see eq 15). All energies are in kilocalories per mole.
v
TABLE 11: SCMTF Iterations for Terminally-Blocked Alanine" iteration no. p i n p i n Xmin 1 2 3 4 5 6 iteration no. 1 2 3 4 5 6 iteration no. 1 2 3 4 5 6
-73.5 -79.7 -79.7 -79.7 -79.1 -79.7
79.2 16.3 16.3 16.3 76.3 76.3
df
df
-74.3 -87.0 -87.6 -87.1 -87.2 -87.6
72.2 106.6 104.1 105.0 106.3 104.9
61.7 60.6 60.6 60.6 60.6 60.6 XC
-39.4 59.4 59.7 61.1 60.5 60.0
Pax
Pax
Xma=
-74.6 -85.3 -85.9 -85.4 -85.4 -85.4
83.1 99.4 97.2 99.4 99.4 99.4
61.1 61.1 61.1 61.1 61.1 61.1
T was taken as 293 K with Mc = 10. The dihedral angles 8, = 4, O2 = +, and 8, = x (in degrees) were obtained by using different convergence criteria, see eqs 19-21. Local minimization from any of these sets of dihedral angles gives the global minimum conformation: 4 = -84O, = 79O, x = 61O. (These values differ slightly from those in ref 44 because of subsequent changes in the ECEPP/2 parameters.)
+
I1 contains the related values of Omax, Omin, and Bc. It can be seen that Omax, t 9 ~ and , Be differ quite significantly. Nevertheless, each of the convergence criteria applied led to the unique global minimum even a t T = 293 K. Computation time was negligible on the IBM 3090 computer. Met-Enkephalin. For Met-enkephalin, the SCMTF iterations were started from 10 different points of the conformational space. These regions are described in terms of the notation of Zimmermann et (1) the energy-minimized a-helical conformation (starting energy, +10.6 kcal/mol), conformation, AAAAA; (2) arbitrarily-selected distorted a-helical conformation originating from a local minimization of conformation AB*AAA (starting energy, 18.6 kcal/mol); (3) seven randomly-chosen conformations with starting energies in the range from +10 to +lo0 kcal/mol; (4) the global-minimum conformation7J6 (FDCIBE) with energy of -13.6 kcal/mol. The first macro-iteration of the SCMTF method was carried out by using a temperature-dependent density for T = 293 K. For all starting points used, the structures that were found in this macroiteration can be described in the notation of Zimmermann et al.4s as FXYEE, where XY = CD*, DC*, DD*, CC*, with the energy in the range of -1 to -4 kcal/mol. When this result is compared with the conformational code of the global minimum (FDC*BE), it can be seen that, whereas the global minimum was not attained, nevertheless these conformations are already relatively close to the global-minimum one. Further, about half of the dihedral angles are very close (at
+
(45) Zimmermann, S. S.; Pottle, M. S.; Nemethy, G.; Scheraga, H. A. Macromolecules 1977, IO, 1.
4616 The Journal of Physical Chemistry, Vol. 96, No. 11, I992
Olszewski et al.
TABLE 111: SCMTF Macroiteration for Met-Enkephalin Starting from the Global-Minimum Conformationo iteration oman omin p i n omin 6Yn O r no. 5 sYn
ein
1 2 3 4 5 6 7 8 9
IO
-82.0 -82.0 -80.8 -80.8 -80.8 -80.3 -79.7 -77.5 -8 1.4 -80.8
149.6 148.5 146.8 144.5 144.5 142.8 147.3 144.5 146.2 146.8
-143.9 -141.7 -137.2 -141.1 -143.9 -155.2 -151.3 -165.9 -165.9 -166.5
82.5 78.0 80.3 88.7 76.3 82.5 65.6 69.6 62.8 56.6
89.3 89.9 89.9 110.7 96.6 93.8 87.0 119.7 133.8 114.1
-73.5 -70.7 -66.8 -51.0 -58.3 -61.7 -52.7 -36.3 -28.5 -22.3
-138.9 -144.5 -134.4 -152.4 -154.6 -154.6 -155.2 -152.4 -153.0 -160.3
23.4 23.9 25.1 153.0 20.6 -170.4 166.5 162.0 157.5 160.3
-165.4 -161.6 -163.6 -163.1 -164.2 -156.9 -165.9 -167.6 -165.4 -164.8
omin 10
conformational codeb
164.2 164.8 164.8 164.8 165.9 164.8 166.5 165.4 165.9 165.4
FDC*DE FDC*DE FDC*DE FDD*EE FDC*BE FDC'EE FDC*EE FDD*EE FDD*EE FDC'EE
TABLE I V SCMTF Macroiterstion for Met-Enkephalin Starting from One of the Conformations Obtained at 293 K' iteration emin Omin 6;ln Oy'n plln 6y gmin emin gmin no. I 5 7 8 6y 10 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1SC 16
litd
-78.6 71.8 -78.6 -78.6 -78.6 -78.6 -78.6 -83.7 -83.7 -83.7 -83.7 -83.7 -83.7 -83.7 -83.7 -83.7 -86
156.3 156.3 142.8 142.8 142.8 142.8 158.0 158.0 158.0 158.0 156.3 154.1 154.1 154.1 154.1 154.1 156
-146.2 -154.6 -157.5 -157.5 -157.5 -157.5 -148.5 -154.6 -154.6 -154.6 -154.6 -154.6 -154.6 -154.6 -154.6 -154.6 -156
169.3 -62.0 156.3 156.3 156.3 -74.6 -74.6 79.2 79.2 79.2 79.2 79.2 79.2 79.2 85.4 85.4
101.7 101.7 79.2 79.2 79.2 79.2 79.2 84.2 84.2 84.2 84.2 84.2 84.2 84.2 84.2 84.2
-101.1 -135.5 -67.3 -67.3 -67.3 -67.3 -67.3 -74.6 -74.6 -74.6 -74.6 -74.6 -74.6 -74.6 -74.6 -74.6
83
84
-74
-156.9 -157.5 -157.5 -147.9 -145.1 -154.6 -146.8 -137.7 -137.7 -137.7 -137.7 -137.7 -137.7 -137.7 -137.7 -137.7 -137
160.3 158.0 21.1 21.1 21.1 21.1 21.1 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 19
-165.4 -164.8 -163.1 -163.1 -163.1 -163.1 -163.1 -165.9 -165.9 -160.8 -160.8 -158.0 -158.0 -164.8 -164.8 -164.8
167.0 167.0 166.5 166.5 166.5 166.5 166.5 166.5 167.0 167.0 164.8 164.8 164.8 164.8 164.8 164.8
-164
160
conformational codeb FEC*EE F*GC*EE FEC*DE FEC*DE FEC*DE FGC*DE FGC*DE FDC*BE FDC*BE FDC*BE FDC*BE FDC'BE FDC*BE FDC*BE FDC*BE FDC'BE FDC*BE
Reference 45. 'The global minimum of the energy a For these iterations, T = 0 K was assumed. See footnote a of Table 111 for description of 6"s. was attained in iteration no. 15. dValues from ref 8 (slight deviations from these values are due to round-off of the values of w's and x's from ref 8, and to the change of the type of the carboxyl hydrogenatom).
most So difference) to those corresponding to the global minimum. This may mean that the global minimum arises when many of the dihedral angles are determined by favorable short-range interactions, whereas some other dihedral angles change in order to satisfy long-range interactions. Even when the global minimum was chosen as the starting conformation,'J6 convergence was not achieved (see Table 111), but Bmin oscillated within loo at the border of the FDD*EE and FDC*EE regions instead of returning to the global minimum. We conclude that, at T = 293 K, the global-minimum potential well is too MITOW to produce any stable state of enkephalin, but instead the molecule chooses a rather very wide region of the conformational space with the flexible glycines changing their conformations over a wide range. When the energy of the structure obtained in the last iteration was minimized, the conformation having a type I' 8-bend with Gly2 and Gly3 at the corners of the bend was obtained. The conformational code for this minimum is FDB*EE and its energy is -3.04 kcal/mol. The diffusion equation method also indicates the presence of a very narrow global-minimum potential well and some much broader wells connected to changes of the dihedral angles of the two glycines.46 In order to reach the global minimum of the energy, it was necessary to lower T and (to improve the convergence) to mix the density distributions from a few preceding iterations. When T was taken as 0 K, and any of the above conformations with (46) Kostrowicki, J.; Scheraga, H . A. J . Phys. Chem., submitted.
energies in the range of -1 to -4 kcal/mol was used initially, the procedure converged to the global-minimum structure, Table IV. Usually up to 10 macroiterations are enough to locate the global minimum exactly. On the other hand, the global minimum could be obtained much more easily, even at higher temperatures (e.g., 293 K), when we rescaled (increased by a factor of 100) the moments of inertia Zi in order to assure that the probability densities of the dihedral angles are sufficiently narrow. One macroiteration for ENK required about 5-10 min on one processor of the IBM 3090 computer depending on the rate of convergence. This may be compared with up to 10 h on six processors of the same machine that were needed for the MCM8 or EDMC'O methods with the side-chain dihedral angles varied (24 variables). Acknowledgment. This work was supported by research grants from the National Institutes of Health (GM-14317), from the National Science Foundation (DMB84-01811 and DMB 9015815), and from the National Research Council of Poland. The computations were carried out a t the Cornell Supercomputer Facility, a resource of the Cornell Center for Theory and Simulation in Science and Engineering, which receives major funding from the National Science Foundation and the IBM Corporation, with additional support from New York State and members of its Corporate Research Institute. We thank Professor B. Jeziorski for reading and commenting on the manuscript and Dr. P. Cieplak for many helpful discussions. Registry No. NANMA, 19701-83-8; ENK, 58569-55-4.
