J . Phys. Chem. 1992, 96,1442-1449
7442
Application of the Diffusion Equation Method for Global Optimization to OMgopeptides Jaroslaw Kostrowickit and Harold A. Scheraga* Baker Laboratory of Chemistry, Cornel1 University, Ithaca, New York 14853-1 301 (Received: April 3, 1992)
The diffusion equation method (DEM) for global optimization, previously applied to simple mathematical functions and to clusters of argon atoms, is modified here for a conformational analysis of two peptides, terminally-blocked alanine and the pentapeptide Met-enkephalin. In this modification,the diffusionequation is solved in the space of the Cartesian coordinates of all atoms of the molecule, and the solution is then examined on the manifold that corresponds to fixed bond lengths and bond angles [a constraint that is characteristic of the ECEPP (Empirical Conformational Energy Program for Peptides) potential function]. The DEM involves a deformation of the potential surface until a "time" tor at which only one minimum remains, and then a time-reversing procedure to f = 0, where the original surface is attained. The modified method was able to locate the global minimum for terminally-blocked alanine, and regions near the global minimum for Met-enkephalin, in a very small amount of computer time [ 0 because such penalty functions evolve differently. For a penalty function that is finite for zero distance between atoms and tends to infinity as this distance increases, the molecule will collapse; Le. as such a penalty function evolves in time, its minimum shifts to smaller interatomic distances at larger times and drives the whole molecule to a collapsed configuration. This is demonstrated in Appendix A for the quartic penalty function, k(?-ro2)2,for which it is easiest to solve the diffusion equation. For the penalty function in expression 4, a collapsed structure will also arise but generally at a different time. If the penalty function does not tend to infinity, the molecule can even expand; this would happen for a penalty function in the form of a linear combination of appropriate Gaussians (see eq 12 of reference 5 for details). Clearly, there would be an arbitrariness in the choice of a penalty function when solving the diffusion equation with the addition of such a penalty function to E(rl,-,r,). Thus, the result would depend on the behavior of the penalty function at large deviations from the assumed geometry, which is nonphysical because the penalty function is physically defined and harmonic only for small deviations. Therefore, it is not justified to solve the diffusion equation in the whole space R3“with the boundary condition 4. The above arguments suggest that a much more justified approach would be to fix the bond lengths and bond angles at acceptable values, use torsional angles as the independent variables, and solve the diffusion equation for an energy E that is a function only of torsional angles. This approach, however, would involve an enormous mathematical difficulty. In the earlier application of the DEM to argon clusters? the diffusion equation in Cartesian coordinates was solved analytically by making use of a FourierPoisson integral. However, evaluation of a Fourier-Poisson integral for a painvise interatomic potential that is expressed as a function of torsional angles cannot be carried out analytically. Consequently, we have elected to solve the diffusion equation in the space of Cartesian coordinates of all the atoms and then consider the solution only for those geometries that satisfy the imposed conditions of fixed bond lengths and bond angles (see section 2). This means that, among all possible arbitrary evolution patterns of the manifold (corresponding to minimal values of the penalty function), we have chosen to use one in which this manifold is time-independent, i.e., in which the values of the bond lengths and bond angles are time-independent. 2. General Approach As was already pointed out in the previous section, we decided to consider the energy as a function of the Cartesian coordinates of all the atoms constituting the molecule. This means that eq 1 with the boundary condition (2) was solved analytically. Then, all the formulas for the energy E(rI,-,r,,),for the given force field, were replaced by the corresponding formulas for the solutions of the diffusion equation, with time t as a parameter. Routine procedures to calculate Cartesian coordinates from dihedral angles were then adopted to enable us to consider the solution only for those geometries satisfying the conditions of fued bond lengths and bond angles. Hence, let @(B,,-,B,,,) denote the vector of Cartesian coordinates of all atoms for the m torsional angles Bl,-.,Bm. @ defines the manifold of fixed bond lengths and bond angles. The function that is minimized locally for a given time t in our approach is F(elr...,emrt) = ~[e(e,,...,e,,,),t]
repeated until t = 0 is reached. Of course, “time-rwersings does not mean that diffusion is a reversible process, i.e. that we solve the diffusion equation backward. We only take the solution forward for the energy function for several time values that form a decreasing sequence. In our experience thus far with the DEM, we have never encountered the problem of branching, Le., where the timereversing trajectory splits into two paths. However, it is possible for a trajectory to end before reaching t = 0, and this problem was discussed in section 2 of ref 5.
3. Modification of the Potential Function We began by using the ECEPP (Empirical Conformational Energy Program for Peptides) potential function”’ to describe the energy of a polypeptide. Bond lengths and bond angles are fixed in ECEPP, and the independent variables are torsional angles. The ECEPP potential can be expressed as A,
-rL2+ - +rhQ- + -Cij r$
Dij rij
+40,
(7)
where ri, is the distance between atoms i and j , and E,, is a torsional energy discussed below. The terms in inverse powers of r are infinite at zero distance. This may introduce an undesirable influence on the competition between minima that are close to regions of high energy. Therefore, to solve the diffusionequation analytically, and to avoid this undesirable influence as r 0, we truncated the interatomic interaction potentials by keeping them constant (and finite) for values of rij less than a certain yi, that depended on the types of interacting atoms i and j (see below). Since it is simpler to obtain an analytical solution of the diffusion equation for truncated Gaussians rather than for the truncated functions r-I2, r4, and f l 0 , we used a linear combination of truncated Gaussians instead of the 12-6 and 12-10 potentials. Thus, the nonbonded and hydrogen-bonded potentials (12-6 and 12- 10, respectively) were expressed as
-
where (9)
and L is the number of Gaussians used. The parameters ui and 6: were calculated separately for each type of (id9 pair of atoms from the parameters ak and bk of the best fit sum-of-Gaussian approximationsfor the r-I2- r4 or (I2 - r-1° functions (see Table I). By this means, we took advantage of the fact that the functions -1,2 - j3r4 or a V 2- b V o can easily be rescaled to the functions r-12 - rd or r-I2 - r-Io, respectively (see footnote 20 of ref 5 for details). yij was chosen as yij = adoijwhere d o , is the distance between atoms i and j for which the potential is zero. The parameter ct was set to 0.75 for 12-6 nonbonded and to 0.77 for 12-10 hydrogen-bond interactions. We decided to make the truncation radius y i . dependent on the radii of interacting atoms to avoid large difkrences in truncated energy values and to make sure that the distance for which interacting atoms have a minimal energy value is larger than the truncation radius. The values 0.75 and 0.77 arise from the requirement to have small truncated energy and the smallest possible truncation radius. The electrostatic interaction between atoms i and j was represented by a truncated Coulombic potential
(6)
We apply the timereversing procedure to the function F (see refs 3-5 for details). This means that we start with a large time value for which we expect F to have only one minimum. We locate this minimum and then we decrease time slightly and minimize F again, starting from the minimum found previously. This is
Bij
&](rij)
= qi4jHy&rij)
where llr
for
r>y
l/y
for
rsy
(10)
Kostrowicki and Scheraga
7444 The Journal of Physical Chemistry, Vol. 96, No. 18, 1992
The parameter yijwas taken to be the same as the corresponding NB parameter (or as the corresponding HB parameter). The term E,,, in eq 7 has the form E,,, = uk(i+ COS P k w (12)
Y
c
(bonds) k
where B is a torsional angle and p is the periodicity. This function is not easily convertible to a suitable form for an analytical solution of the diffusion equation. The cosine of a torsional angle can be expressed as cos e k ak + bkrz(1,4)k (13) where r(l,4)k is the 1,4 distance for the kth torsional angle and ak,& are constants calculated for fned bond lengths and bond angles. Generally, the cosine of a multiple torsional angle, cos can be represented as a Tschebyscheff polynomial in cos 0, likewise as a polynomial in r(1,4?, and can easily be transformed according to the diffusion equation. This suffers, however, from the major disadvantage that eq 13 is valid only on the manifold of fmed bond lengths and bond angles, and the right-hand side of this equation very rapidly tends to --03 off this manifold because b k in eq 13 is always negative (on the manifold, the distance r(1,4)is lowest for B = 0 and increases as B increases toward r). During the time evolution, this off-the-manifold region of very low function values has an undesirable influence on the situation on the manifold. It appears that, instead of using eq 13, one should use a form that reproduces E,, correctly at physically-meaningful distances between bonded atoms and decays rapidly to zero for other distances. Such a form might involve the exponential boundary condition of eq 5 to transform EIor. We are currently exploring the use of such an exponential boundary condition. In the present preliminary application of the DEM to oligopeptide8,we have decided to circumvent this problem temporarily by excluding the term Em, i.e., by setting the rotation barriers uk to zero for the side-chain torsional angles x (these barriers are small, and we are interested here primarily in identifying the region with the best backbone conformation). For the backbone torsional angles 4 and the barrier is actually zero in ECEPP. For the torsional angle o about the peptide bond, the bamer is large, but fortunately these torsional angles rarely depart substantially from 180O; hence, we have fmed all a’s at this value. The diffusion equation was solved analytically for truncated Gaussians and a truncated l / r function in the space of Cartesian coordinates of all interacting atoms. This means that the functional forms of the expressions for pairwise interaction energies were replaced by the formulas in Appendix B. 4. Estimation of the Starting Time t o The initial time to for starting the time-reversing procedure should be sufficiently large so that only one minimum remains in the deformed potential surface. The solution, 6(rl,--,rn,r), can be expressed as 6(rI,-,rn,t) =
v),
+,
Figure 1. Illustration of the most distant conformations of two chains. The open circles, numbered 1, 2, -, are the atoms of the first conformation, with their cOOrdinates being rl, r2, -. The open circles,numbered l‘, 2’, -, are the atoms of the second conformation, with their coordinata being pl, p2, -. In each of these two conformations the bond lengths are not nectssarily equal. The filled circles are the positions of the atoms in the fully-distorted conformations, directed opposite to each other and with all bond lengths set equal to the largest one in the original chain. The half-filled first three atoms have the same fixed positions in all conformations.
are any two geometries satisfying the constraints of fixed bond lengths and bond angles and situated in the manner described (hence, lying on the manifold) in the coordinate system. The distance 6 between these two geometries cannot be larger than the value adopted by an arrangement of two fully-distorted chains directed opposite to each other (with the first three atoms of each coinciding) and with all bond angles taken as 180° (except for the first ones in both chains which have the original value and the second one in the first chain which is taken as O O ) ; see Figure 1. Let d be the largest bond length in the chain (it is justified to consider only the largest when an upper bound is being established). Then an upper bound on 62 of eq 15 is given by
6’ < (24’ One of the two minima of E that lie closer to each other than the width K = 2 ~ of~the1exponential ~ part of the integrand in eq 14 will usually disappear before time t (seeAppendix C for details). Hence, when the width K is larger than the largest dimension of the manifold on which we search for the minimum (the manifold corresponding to fixed bond lengths and bond angles), we can expect that only one minimum will remain on the manifold. We may estimate the largest dimension of the manifold as follows. In our calculations, the first atom is situated at the center of the coordinate system, and the second atom is on the x-axis. The largest dimension of the available part of the manifold (that part in which the first and second atoms are fixed) is the largest value of the expression 6(r1,-r,, Pl,.*.#n) = (Zlh - pk112)1’2 (15) where eq 15 is an RMS-like quantity, and (rl,-.,rn) and ( p l , * ~ , )
2 + (4d)2 + + (2nbd)2 = -d2nb(nb + 1)(2nb + 1) 3 (16)
+
where n b 2 is the number of bonds in the chain. Equation 16 was obtained by assuming that each bond has a length d, and summing over all pairs of corresponding atoms (eq 15) [in Figure 1, r1 = pI, rz = PZ. r3 = ~ 3 Ilr4 , -Pd2= (2d2, llrs - ~ 5 1 1 ’ = (44’, etc. for the chains with filled circles]. The condition for obtaining to requires that 6’ < K2
This inequality is satisfied for 1 to
> gdlnb(nb
4to
+ I)(& + 1)
(17)
(18)
The inequality (18) should be considered only as a rough estimate because the proof in Appendix C is based on a particular form
The Journal of Physical Chemistry, Vol. 96, No. 18, 1992 7445
Global Optimization of Oligopeptides
TABLE I: Gaussian Parameters of the Best-Fit Functions to Two Standard Potentialsa a' b' a2 b2
-1.63 1 1.353 2008.734 8.468
-1.128 2.148 2312.359 9.773
Third-Component Parameters
a3 b3
-1 .o 1o4
-1 .o 1o4
OAs shown in footnote 20 of ref 5, these parameters can be rescaled for any type of i j pair of atoms.
of a minimized function and the estimate of the distance 6 given by eq 16 neglects the side-chain atoms. Having estimated to, it is necessary that the steps in the time-reversing procedure be small for small values of t because the solution of the diffusion equation changes more rapidly as t approaches zero. Hence, we adopted the following protocol for selecting the time values in the time-reversing procedureS tq = to[1
-
];
-I 80
-90
0
90
I80
4 (deg)
Figure 2. Deformed potential surface for terminally-blocked alanine at t
= 10
A2. The position of the unique minimum is indicated by M.
2
where q is the step number and Q is the total number of steps.
5. Some Details of tbe Algorithm After to was estimated, the formulas in Appendix B were used instead of the original (ECEPP) expressions for the energy for various values oft < to. The starting conformation corresponded to the unique minimum at to. At each time step, the function F(O,,-,O,,t)of eq 6 was minimized locally, starting with the values of the torsional angles obtained in the previous time step. (6 of eq 6 is the sum of 9 and % of Appendix B.) In this minimization, the Cartesian coordinates of all atoms (rl,--,rn) were computed from the corresponding set of torsional angles, assuming that the bond lengths and bond angles are fixed, and inserted into the analytical expression for the solution 6(rl,-,rn,f) of the diffusion equation in Cartesian space. In this way, it is guaranteed that the solutions in Cartesian coordinates correspond to the manifold on which the bond lengths and bond angles are fixed. These Cartesian coordinates were also used to calculate the derivatives required for the minimization (although they are derivatives with respect to torsional angles, they are expressed in terms of Cartesian coordinates). After a sufficient number of time steps to reach t = 0, the values of (O,,--,Om) and F(O1,-,8,,O) are expected to correspond to the global minimum on the manifold. 6. NumericalResdts
Two small peptides were chosen to test the methodterminally-blocked alanine and the pentapeptide Met-enkephalin, H,N-Tyr-Gly-Gly-PheMet-COOH. The first molecule, consisting of 22 atoms, with only the torsional angles 4 and $ as variables, was used as a test to check whether the solution of the diffusion equation in a high-dimensional space of Cartesian coordinates of all of the atoms is smoothed properly on a relatively low-dimensional manifold. The second molecule (consisting of 75 atoms) was selected as a sufficiently complex system whose global minimum was known from previous extensive computation^,'^-'^ and it can be assumed that all of the deepest minima had been explored previously. Only 19 vEriables were chosen for this problem (4, $, and x's of each residue). Since the number of local minima for this molecule is huge12 (more than lo"), this is a more stringent test of the DEM as a global-optimization algorithm. Terminally-Blocked Alanine. The nonbonded and hydrogenbonded functions of ECEPP were represented by a linear combination of two truncated Gaussians [eqs 8 and 9 and Table I (L = 2)] and the electrostatics by a truncated l/r function (eqs 10 and 11). Since only the torsional angles 4 and $ were varied, with the three CH3 groups being fixed at their global minimum pos-
-180
- 90
0
90
I80
4 (deg)
Figure 3. Deformed potential surface for terminally-blocked alanine at
t = 0.295 A2. The unique minimum of the t = 10 A2map of Figure 2 has moved to the intermediate position I but (even though other minima may appear) there is no problem in identifying I with the position of the global minimum of the original ECEPP function, Le., at t = 0. The contours in the lower right-hand corner correspond to a maximum.
itions,l* the omission of E,,, did not affect the computations. For the time to = 10 A2, the (4,$) map exhibits only one minimum, in the upper left-hand corner (Figure 2). There is a maximum at 4 1 5 O , $ - 1 5 O and saddle points at 4 loo, $ 165O and at 4 165O, $ -10'. An example of a (&+) map at an intermediate time (t = 0.295 A2) during the time reversal is shown in Figure 3, from which it can be seen that, even though other minima reappear, there is no trouble identif ing the shifted position of the unique minimum of the t = 10 map. As t 0, this minimum follows the trajectory depicted in Figure 4, and ends up as the global minimum on the original (&+) map for t = 0 (Figure 5). Met-enkephalin. Substitution of d = 1.53 A and nb = 15 in eq 18 gives the following rough estimate for to for Met-enkephalin
-- - -
-
-
w2
-
to > 2.9 x io3 A2
(20)
Thus, to be overly cautious, we took to as 104 A2. We started with the nonbonded and hydrogen-bonded functions of ECEPP, r e p resented by two truncated Gaussians, and with a truncated ECEPP 1 /r function for electrostatics (as for terminally-blocked alanine). At to = lo4 A2,about 100 randomly-chosen conformations were selected for local minimization on the deformed potential surface. Since all minimizations converged to the same minimum, we may conclude that the value of to = lo4 A* was sufficiently large to smooth the surface to the required extent.
Kostrowicki and Scheraga
7446 The Journal of Physical Chemistry, Vol. 96, No. 18, 1992
TABLE II: Computed Conformations of Met-enkephalin torsional angles (deg) Tyr-1 Gly-2 Gly-3 Phe-4 Met-5 Tyr-1 Tyr-1 Phe-4 Met-5 Met-5 confor- 4 4 4 4 4 x' x3 x' x' x3 EOsb (1/2)CihWi A$b AQ &b mation $ (t $ $ $ x2 x2 x2 x4 (kcal/mol) (kcal/mol) (kcal/mol) (kcal/mol) (kcal/mol) -68 -168 -170 85 -93 -94 -138 173 -74 A -83 -1.60 1.4 -18.2 -19.8 -18.2 ,
B
c
-104 -166 62 -160 71
151 -87 155 -86 155
72 81 -77 80 -74
156 -134 169 -157 159
152 -102 55 -117 12
93 64 69 76 79
116 65
83 -55 102 62 91
-172 -62 -168 49 168
-55 -175 173 -176 -51
-8.06
2.2
-14.0
-22.1
-14.1
. -13.50
2.6
-11.6
-25.1
-1 1.7
a All energies correspond to a two-Gaussian representation of the nonbonded and hydrogen-bond interaction and differ slightly from the original ECEPP ~ a l u e s . ' ~ -bThis ' ~ is the same as F(8O) of eq 26.
L
3
- 90
0
90
d (deq 1 Figure 4. Trajectory of the global minimum for terminally-blocked alanine. Point M corresponds to the unique minimum at r = 10 A*, and point C corresponds to the global minimum at r = 0.
on argon cluster^,^ we circumvented this problem by adding to the two-Gaussian representation of the Lennard-Jones potential a third, very wide Gaussian of relatively small magnitude; this forced the potential to evolve into a purely attractive one. This additional Gaussian was found to be important for obtaining low-energy minima for argon clusters. Consequently, in a similar manner, we applied three truncated Guassians (with parameters given in Table I; L = 3 in eq 8) as a representation of the nonbonded and hydrogen-bonded functions (with the electrostatic function remaining unchanged). With three Gaussians, the interatomic potentials were attractive for to = lo4 A2. This suggested the possibility that more than one minimum might exist at that time. An extensive minimization study, however, revealed that only two different structures of the backbone are present at to = lo4 A2 for the three-Gaussian form of the potential. Unlike the conformation with the two-Gaussian potential, the two structures were rather compact. The trajectories of both minima led to (the same) structure B (Table I1 and Figure 6) at t = 0 which is practically the same as the reported global minimum (structure C of Table 11). The backbone of structure B exhibits the &turn that is characteristic of the lowest-energy structure C. The differences in the side-chain conformations are due to the neglect of E,,,. 7. Statistical Weights and the DEM
It is well-known that the free energy (or statistical weight) rather than the energy determines the most stable conformation. Consequently, in order to interpret the results in Table 11, it is necessary to take account of the vibration of the molecule in different conformational states. Although we have analyzed the conformational states of a chain molecule with fmed bond lengths and bond angles in the present study, there are still many degrees of freedom for motion arising from changes in the torsioinal angles.16 Even at T = 0 K, the energy of such a chain is not the value at the minimum of the potential; it is higher by the zero-point librational energy per0 = (1 /2)ChWi (21) I
-I 80 -I 80
- 90
0
90
I80
4 (deg) Figure 5. ECEPP (+*$) map of terminally-blocked alanine at r = 0. Point C is the position of the global minimum.
The time-reversing procedure yielded backbone-plus-side chain conformation A of Table 11. Its energy was relatively high (-1.6 kcal/mol) compared to the value of -13.50 kcal/mol for the global-minimum conformation (we minimized the Gaussian form of the potential, starting from the global-minimum conformation found in ref 12; it converged to conformation C of Table 11, with an energy of -13.50 kcal/mol). Failure to find the global minimum when using the two-component truncated Gaussian representation is probably due to the fact that the potential does not evolve into a purely attractive one at large t. Because of its large positive value at small distances, it becomes repulsive at large t for all distances attainable for pairs of atoms in the molecule. This is the reason that elongated shapes are preferred at large t (not shown here). In previouS computations
where the oi)sare the normal-mode frequencies. The squares of these frequencies, a?,can be obtained, in a harmonic approximation, as the eigen~alues~~J* of the matrix r13, where 3 is the matrix of second derivatives at a minimum (eol, eo,) -0,
and r is the matrix for the kinetic energy T of the internal librational motion of the molecule, expressed in terms of time derivatives of the torsional angles. I' is defined by the formula ij
To separate the internal motions from the overall translational and rotational motion, the Eckart conditions were applied.'*J9 The calculated zero-point librational energies for minima A, B, and C of Met-enkephalin are presented in Table 11, from which it can be seen that this correction at T = 0 K is not the same for each conformation.
The Journal of Physical Chemistry, Vol. 96, No. 18, 1992 7447
Global Optimization of Oligopeptides
d
Figure 6. Stereoview of conformation B of Met-enkephalin (Table 11). The thin line represents a hydrogen bond.
At temperatures above absolute zero, the molecule can exist in excited;ibrational states. The quantum-mechanical partition function for the vibration can be expressed as16a e
b
= n [ 2 sinh (@hai/2)]-I I
(24)
The Helmholtz vibrational free energy A$b
= -(l/@) In
e
(25)
b
at T = 298 K for conformations A, B, and C (where @ l/kBT) is presented in Table 11. The total free energy of the molecule is the sum of the energy at the minimum and the vibrational free energy correction, viz.
AQ = F(8')
+ A$b
(26)
The values of AQfor conformations A, B, and C (given in Table 11) are closer to each other than are the corresponding minimum energies. It can be seen that A$,makes a significant contribution to the result obtained by the DEM. This behavior can be understood by considering the following discussion of the classical limit of the partition function,20viz.
ZC' = (const)le-fiH@*9) dp dq
(27)
where H(p,q) is the Hamiltonian of the molecule and p,q are the generalized momenta and coordinates, respectively. When the mnformation of the molecule is expressed in terms of the Cartesian coordinates of the atoms, ri, and its motion by the corresponding momenta, pi = mii,integration over the momenta gives a constant, independent of the positions ri, and the partition function takes the form = (const)le-fi~(rI..".4) drl
...dr,
(28)
When some geometrical constraints are imposed on the shape of the molecule, other generalized coordinates, q, rather than Cartesian coordinates, are more convenient.16 If bond lengths and bond angles are fixed (as in this work), the most convenient coordinates are torsional angles, 8. Integration over the conjugate momenta and the external coordinates now yields16
ct,= (const)
[det r(8)] 1/2e-fiF(e) d8
(29)
The classical formulas, eqs 28 and 29, are valid only for sufficiently high temperature. To determine the difference between the classical and quantum approaches for Met-enkephalin, we calculated the vibrational free
energy also in the classical approximation. By inserting the harmonic form16aof F(0) F(8) = (1 /2)CSij(Oi - Ooi)(Oj ij
- eoj)
(30)
into eq 29, and replacing r(O)by the value at the minimum, r(Oo), the following formula is obtained for the classical vibrational Helmholtz free energy?
A$, = (1/2/3) In [det(@2h21+13)] (31) The calculated values of A$: for conformations A, B, and C are presented in Table 11. Since the quantum and classical vibrational free energies are in very good agreement, we base our further discussion on the classical treatment only. There is, however, one point that needs further discussion, viz., the existence of detr(8) in the integral in eq 29. It was argued previously16bthat the hard vibrational modes (e.g. bond stretching) vary with changes of conformation in such a way that this determinant is compensated by changes in the frequencies of the hard modes and, therefore, that a more accurate treatment for the free energy is given by the integral without this factor. This conclusion was arrived at, however, for a flexible model with very large force constants for bond stretching and bond-angle bending. This model is more physically realistic than one in which the bond lengths and bond angles are fixed exactly. However, we are concerned here primarily with the mathematical relation between the partition function for a well-defined system (i.e. one with fmed geometry) and,the time evolution of the potential according to the diffusion equation. When the geometry of the molecules is exactly rigid, then eq 29 is the correct expression for the partition function, as pointed out in ref 16b. The factor [detr(0)]1/2 has a clear geometrical interpretation; i.e. [detr(0)]1/2d0 is a measure on the manifold of those conformations for which the bond length and bond angle constraints defined by the function @(e)in section 2 are satisfied. Hence, when the function e-fiE(rl*-sdor e ~ F ( e l *is- *applied e~ as the boundary condition for diffusion in the whole space of Cartesian coordinates or on the manifold, respectively, the integrals in eqs 28 and 29 are analogous to the initial balance of heat or mass. Indeed if one considers the function e-@E(rl*".JJ as the spatial concentration distribution of a substance in a solution, then the integral (28) will be the total amount of that substance. This amount is consewed throughout the physical process of diffusion of the substance in the solution. Mathematically, this means that one can integrate the concentration distribution for any time to obtain the same value of the integral as for the initial distribution. Hence, the partition function calculated from the solution of the diffusion equation
7448 The Journal of Physical Chemistry, Vol. 96, No. 18, 1992 will be the same for any time, provided that the boundary condition
is applied in the exponential form. When the integration in eqs 28 or 29 is carried out over only the part of the space or manifold that corresponds to a certain energy minimum, the integral can be interpreted as the statistical weight of that minimum. When the integration is confined to the region around the minimum where the energy can be assumed to be harmonic, then the integrand will be a multidimensional Gaussian. We proved earliers that, during the time evolution of well-separated Gaussians, the one of the largest integral will dominate at large times. Hence, the minimum that survives and is extracted by the DEM is the one with the largest statistical weight (largest probability). It should be stressed that this conclusion is based on several assumptions, the most important being that the classical expression for the statistical weight is valid and that the potential energy is transformed exponentially prior to solving the diffusion equation. Use of the potential energy itself, instead of the exponential form, as the boundary condition to solve the diffusion equation, as is done in this paper, should be understood as the high temperature approximation when e-BE does not differ much from 1 - BE. Another comment is pertinent to the method of minimization with constraints. The above comments show that the results of the DEM can generally be interpreted physically only when diffusion proceeds on the manifold, so that the conservation equation (29) is satisfied. Reference was made to this point in section 1.
8. Computing Time It is of interest to cite the computing time for one processor of an IBM 3090 computer. No computing time is involved in “reaching” to; Le. the value of to is simply estimated by the procedure in section 4. For terminally-blocked alanine, a (&$) map at to, obtained in - 5 min, indicated the presence of only one minimum; the time-reversing procedure involved 100 steps for a total CPU time of -1 min. To test that only one minimum remains at to for Metenkephalin, about 100 energy minimizations were carried out at this time, requiring 10 min. The time-reversing praccdure involved 300 steps, for a total CPU time of 500 s for this molecule. This is the time required for local minimizations with the Gill-Murray method,2I with first and second derivatives calculated analytically with an algorithm developed in our laboratory. For comparison, the MCMI3 and EDMCI5 procedures required approximately 2-3 h on the IBM 3090 computer (with six processors) to locate the global minimum for Met-enkephalin. However, whereas MCM and EDMC use random starting points, the DEM is a deterministic procedure. The computation time would be expected to scale cubically or quartically with the number of atoms. This estimate is based on the observation that the number of time steps required in the timereversing procedure does not seem to depend very significantly on the number of atoms and that the number of iterations in local minimizations does not increase more rapidly than linearly with the number of atoms. The cost of calculating a function and its derivative increasesquadratically with the number of atoms. The total cost arises from the number of time steps, the number of iterations in a local minimization, and the cost to evaluate a function and its derivatives. The product of the above factors should be of the order of n3 or at most n4.
-
9. Conclusions The DEM succeeded in finding the global minimum for terminally-blocked alanine and the region of the global minimum for Met-enkephalin. It appears, however, that the width of the minima and conformations of large energy in the neighborhood of minima contribute substantially to the competition between minima. Hence, when one is interested in finding the deepest minimum regardless of its width and possible high-energy surroundings, one needs to consider carefully modifications of the potential function prior to applying the DEM. The most desirable modification would be the exponential transformations-*which not only emphasizes low-energy regions
Kostrowicki and Scheraga but also allows one to interpret the results in terms of a temperature dependence. This is, however, difficult to apply mathematically, and is currently under investigation in our laboratory. In this paper, other means of excluding high-energy regions were adopted, e.g. truncation of the potential and addition of a third negative-integral Gaussian. They proved to be successful in searching for regions of lowest energy, although difficulties in dealing with the complete, exact potential function were encountered. Acknowledgment. This work was supported by the National Institutes of Health (GM-14312), the National Science Foundation (DMB84-01811 and DMB90-15815), and the National Foundation for Cancer Research. The computationswere carried out on the IBM 3090 supercomputer of the Cornell National Supercomputer Facility (CNSF), a reSOUTCe of the Come1 Theory Center, 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. Appendix A. Solution of the Diffusion Equation for a Quartic
Penalty Function In order to see how the minimum of the penalty function shifts to smaller bond lengths as it evolves in time, consider the penalty function S(ri,rj) E (11s- rj1I2- ro2)2 (A- 1) where ro is either the strain-free bond length or the 1,3-distance; Le. the penalty function becomes zero as the bond length or the 1,3-distancebetween atoms i and j , llri - rjll, approaches ro. The solution of the diffusion equation for this function is
8(ri,rj,t) = [Ilri - rill2 - (ro2- 20t)I2 + 16rO2t- 160r2 (A-2) where eq A-2 arises from the application of the operator T(r) of eq 4 of ref 3 to the function S(ri,rj). The function 8 is of a form analogous to S. It has its lowest value for an interparticle distance equal to (r$ - 20t)1/2. This means that the interparticle distance decreases with time, starting from the value of ro for t = 0 and reaches 0 for t = rO2/20.
B. Solution of the Diffusion Equation Since we have separated the total energy into its components (eq 8 for NB and HB and eq 10 for el), we can solve the diffusion equation for each component separately. Hence, for the NB, HB part, 0 , of the total energy, 6,the analog of eq 1 becomes Appendix
(B-3) The solution of eq B-1 is (B-4) where
[
erfc( y(1 2[t(l
)
y(1 + 4br) + r ++4bt) - r + erfc( 2[r(l + 4bt)]l/2 4bt)]1/2
ie-by2[
2
erf
(=)+
erf
2f‘/2
(-)I
Y+r 21‘12
(B-5)
The Journal of Physical Chemistry, Vol. 96, No. 18, 1992 7449
Global Optimization of Oligopeptides erf (x) is the error function defined by erf ( x ) = ( 2 / d 2 ) X .fie-” du and erfc ( x ) = 1 - erf ( x ) . Analogously, the solution of the diffusion equation, for the electrostatic function
a2 -aI -Sl(x,t) ax2
where a
=
(2(7rt)’lZ)-’
as the boundary condition, is given by Hence, the condition for two minima to merge is -1 -a(;,?) a2 10 a ax2 Then a2 2 - - 1 0 or a 5 2(2r)’l2 4t Registry No. Met-enkephalin, 58569-55-4.
where
(C-7)
References and Notes
Appendix C. Calculation of the Merging Time for Two Infinitely Narrow Minima Although one cannot guarantee that, after a certain time (in the application of the DEM),only one of two close minima of any function will remain, one can at least estimate a rough value of that time by considering a model function that has very deep minima separated by a certain distance. This estimate is based on our observation that the deeper the minimum, the longer the time after which it merges with a neighboring maximum. Thus, less deep minima will disappear more rapidly. As a function with a pair of very deep minima, consider an “energy” w ( x ) that is a sum of two Dirac delta functions separated by the distance a, Le. w(x) = -[6(x)
+ 6(x - a)]
(C- 1)
The solution of the diffusion equation for w is given by
+
n(x,t) = -(2(Tt)1/2)-ye-x3/4r e-(~a)’/49
(c-2)
Sl has a maximum at x = a/2 provided that W(a/2,?) is negative. When W(a/2,t) becomes positive, the maximum at a/2 which separated two minima of Sl becomes the only minimum of this function
-i
a
a ax
= &-2/4l + -e-(x-a)2/41 x-a 2r
2r
(C-3)
(1) (a) Stillinger, F. H.; Weber, T. A. J. Star. Phys. 1988,52, 1429. (b) Head-Gordon, T.; Stillinger, F. H.; Arrecis, J. Proc. Narl. Acad. Sci. U.S.A.
1991,88, 11076. (2) Zakharov, V. V. Eng. Cybern. (Engl. Trans.) 1970,8,637. ( 3 ) Piela, L.; Kostrowicki,J.; Scheraga, H. A. J. Phys. Chem. 1989, 93, 3339. (4) Kostrowicki, J.; Piela, L. J . Optimiz. Theory Appl. 1991, 69, 269. (5) Kostrowicki, J.; Piela, L.; Cherayil, B. J.; Scheraga, H. A. J. Phys. Chem. 1991, 95,4113. (6) Stillingcr, F. H. Phys. Rev. B 1985, 32, 3134. (7) Yuille, A. L.; Poggio, T. A. IEEE Trans. Parr. An. Mach. Intell. 1986, PAMI-8, 15.
(8) (a) Shalloway, D. In Recent Advances in Global Optimization, Floudas, C., Pardalos,P., as.Princeton ; University Res: Princeton, NJ, 1992; pp 433-477. (b) Shalloway, D. Global Optimization, in press. (9) Momany, F. A,; McGuire, R. F.; Burgess, A. W.; Scheraga, H. A. J. Phys. Chem. 1975, 79,2361. (IO) Nhnethy, G.; Pottle, M. S.; Schcraga, H. A. J . Phys. Chem. 1983, 87, 1883. (1 1) Sippl, M.J.; Nbmethy, G.; Scheraga, H. A. J. Phys. Chem. 1984,88, 6231. . .. (12) Li, Z.; Scheraga, H. A. Proc. Narl. Acad. Sci. US.A. 1987,84,6611. (13) Li, Z.; Scheraga,H. A. J. Mol. Srrucr. THEOCHEM 1988,179,333. (141 Purisima, E. 0.; Schcraga. H. A. J . Mol. Biol. 1987, 196, 697. (19 Ripoll, D. R.; ScheragarH. A. J . Protein Chem. 1989, 8, 263. (16) G6, N.; Scheraga, H. A. (a) J . Chem. Phys. 1969, 51, 4751; (b) Macromolecules 1976, 9, 535. (17) Wilson, E. B., Jr.; Decius, J. C.; Cross, P. C. Molecular Vibrations; Dover Publications, Inc.: New York, 1980; pp 54, 273. (18) Noguti, T.; G6,N. J. Phys. Soc. Jpn. 1983, 52, 3283. (19) Eckart. C. Phvs. Rev. 1935.47. 552. (20) Huang; K. Slbrisrical Mechanics; John Wiley h Sons, Inc.: New York, London, Sydney, 1963; pp 157, 213-220. (21) Gill, P. E.; Murray, W. Math. Prog. 1974, 7, 311.
