J. Phys. Chem. 1993,97, 4416-4420
4416
Potential-Scaled Molecular Dynamics and Potential Annealing: Effective Conformational Search Techniques for Biomolecules Hideki TsujishitaJ Ikuo Moriguchi, and Shuichi Hirono’ School of Pharmaceutical Sciences, Kitasato University. Minato- ku, Tokyo 108, Japan Received: December 8, 1992
An effective conformational search method, potential-scaling molecular dynamics, was described for molecular dynamics (MD) simulations using the explicit solvent water. M D simulations in solution generally yield more accurate results than simulations in uucuo, but conformational search methods including the process of raising the temperature (high-temperature MD or simulated annealing) cannot be applied to a system including a water solvent because these methods affect the whole system and cause unrealistic behavior of water. In the potentialscaled MD, the potential energies involved in the desired degrees of freedom (in our cases, those of solute molecules) are partially scaled down, instead of raising the temperature of the whole system. This method was tested by using a simple model compound, N-acetylglycyl-N’-methylamide(1, Figure l), and proved to be capable of searching the conformational space accurately and more efficiently than the normal MD. Using the idea of the potential-scaled MD and simulated annealing, we designed a potential-annealing approach and examined its usefulness. This method could find the lowest energy regions in the conformational space of the model system, starting with the different conformers appearing during the potential-scaled MD simulation. It is suggested that the two methods described here, potential-scaled MD and potential annealing, have potential as powerful tools in conformational analyses of molecules in aqueous solution, especially in predicting the three-dimensional structures of biomolecules. 0
Introduction Two difficulties weigh heavily in conformational analyses and structure predictions of biomolecules. One is the size of the conformational space, and the other is the presence of many local minima on the potential energy surface. Among many procedures proposed for attacking these problems, one of the most powerful methods is molecular dynamics (MD) simulation. In MD simulation, the actual motion of a system at a given temperature is computed by solving Newton’s equation of motion. The kinetic energy of atoms allows the system to surmount energy barriers and the conformational space can be searched spontaneously. MD simulation at high temperature1,*gives the system higher kinetic energy and increases the ability to surmount energy barriers. Consequently,a larger part of the conformational space can be searched at high temperature than at normal temperature. The simulated annealing method3 consists of more complicate processes. At the heart of this method is an analogy with thermodynamics, specifically with the way that liquids freeze or crystallize or metals cool and anneal. The temperature of the system is first increased in order to search an extensive region of the conformational space, and then the temperature is slowly decreased to a sufficiently low level. Generally, this method can find the lower energy state of a system better than an ordinary gradient energy minimization method. The simulated annealing is now successfully used in various fields, especially in building up three-dimensional structures based on X-rayM or NMR7-9 data. However, there is a serious problem when these methods including the process of raising the temperature are applied to a biomolecular system with explicit solvent water molecules. Raising the temperature of the system (typically up to about 1000 K) would affect not only the solute but also the solvent water molecules. Consequently,the pressure of the system would be extremely increased (under constant volume), or the solvent water would turn into vapor (under constant pressure). Therefore, it is virtually impossible to apply the process of raising
’
Present address: New Drug Research Laboratories, Kanebo Ltd., 5-90, Tomobuchi-chou 1-chome, Miyakojima-ku, Osaka 534. Japan.
ti 1
+
Figure 1. Structure of N-acetylglycyl-N’-methylamide. Q and dihedral angles are indicated by round arrows.
the temperature to the solution system directly. However, MD simulation in a biological system is generally more realistic when water molecules are included than in vacuo. When the system is simulated in vacuo, the properties of a t o m near or at the surface of the system can be distorted by the vacuum boundary condition. This boundary condition also distorts the shape of a nonspherical molecule, since it tends to minimize the surface area.I0 A solvent with a high dielectric constant such as water causes a shielding effect on the electric interactions between charges or dipoles in a molecule.Il In simulations in vacuo, the lack of this effect causes electric interactions that are too strong. In addition, intramolecular hydrogen bonding, which plays an important role in protein structure, is affected by the competitive presence of water molecules surrounding a solute in aqueous solution and many spurious intramolecular hydrogen bonds can be formed during MD simulation in vacuo.Io An example of the importance of the presenceof theexplicit solvent water is found in theliterature by Levitt and Sharon.12 In their work the results from the simulations of a small protein, bovine pancreatic trypsin inhibitor (BPTI), with explicit water and in vacuo were compared. The time-averaged structure of BPTI obtained from MD simulation in water was much more like the structure observed in highresolution X-ray studies than that obtained from MD calculation in vacuo. The amplitudes of atomic vibration in solution were smaller, and fewer incorrect hydrogen bonds were formed in the simulation in solution. MD with explicit solvent water is preferred for these reasons, especially if there is little experimental information available. However, it is noted that ordinary methods
0022-3654/93/2091-4416~04.00/00 1993 American Chemical Society
The Journal of Physical Chemistry, Vol. 97, No.17, 1993 4417
Potential-Scaled Molecular Dynamics for conformational searching including the process of raising the temperature cannot be applied to a system with explicit water. Here we report an effective conformational search technique for a biomolecular system using molecular dynamics simulation, which does not involve explicitly raising the temperature.
180 120
60 n
a
Methods
t-!
The population of conformers in a system in equilibrium at a temperature T is governed by the Boltzmann factor:
E
Q,
..P(-&)
= exP(-&)
exP(-&)
(1)
where H represents the Hamiltonian, Kand Yare the kinetic and potential energy, respectively, and kB is the Boltzmann constant. Because kinetic energy is proportional to temperature, the first term on the right side of eq 1 becomes a constant A. When the temperature is increased by a factor p, eq 1 becomes as follows:
This equation indicates that both raising the temperature by p and decreasing the potential energy by p at constant temperature
have the same effect on the Boltzmann factor (1). This means that a larger conformational space of the system can be searched by scaling down its potential energy as well as by raising its temperature. The advantage of scalingdown the potential energy is that the energy of the desired degrees of freedom can be reduced selectivity, while the whole system is affected by raising the temperature. This idea of scaled potential energy was originally found in the work by Mark et aZ.I3 In the estimation of free energy, they used the idea of searching the configurational space quickly and generating the statistical ensemble efficiently. Here we extend this idea as an effectiveconformational search method using MD with explicit solvent water. When scaling down the potential only of the degrees of freedom involved in the solute molecules instead of raising the temperature of the whole system, one can search the conformational space of the solute molecules efficiently without the bad influencesonthe behavior of thesolvent water. Our calculations in this paper used the potential energy function of AMBER version 4.0.14 This function consists of bond stretching, bond angle, dihedral angle, van der Waals, electrostatic, and hydrogen bonding interaction terms. In order to selectively scale down the potential related to the solute molecules, only the bond angle and dihedral potential terms in the potential energy function were reduced by a scaling factor p. Nonbonded interaction potential terms (van der Waals, electrostatic, and hydrogen bonding) play very important roles in the behavior of water molecules. Decreasing these terms, therefore, would cause serious effects on the arrangement of the solvent water molecules. In addition, because all bond lengths were fixed by SHAKEIS during the MD runs (see below), scaling the bond stretching term would be meaningless. The dihedral potential term is one of the major components of the potential energy barriers. Decreasing bond angle terms causes a molecule to be more flexible and would allow the system to pass through alternative paths while surmounting the energy barriers. Moreover, because the internal degrees of freedom in the TIP3P water model16 used here included the bond stretching interactions only, scaling down the bond angle and dihedral terms produces no effects on the internal motion of the water molecules. Therefore, reduction of these two terms is sufficient to make the conformational changes in the soluteeasier and would not cause any unfavorableinfluences on the solvent. Using the concept of this potential-scaled MD, a method like simulated annealing can be also designed. The system is first simulated long enough under the condition of low-scaled potential
0
-60 -120 -180 -180
-120
-60
0
60
120
180
4) (degree) Figure 2. Potential energy surface of 1. Darker parts represent lower energy regions. Each label (A-E) indicates an energy minimum. energy. This corresponds to the high-temperature state in simulated annealing. Then the potential of the system is gradually increased, equivalent to the slow cooling process in simulated annealing. This is the technique of annealing, and it is essential for ensuring that a low-energy state will be achieved. The final structure obtained is expected to have a sufficiently low energy level. If the parameters are chosen appropriately, this potentialannealing method can be expected to find the structure in the global minimum automatically. The effectiveness of the potential-scaled MD was tested by using a small compound, N-acetylglycyl-Ncmethylamide(1, Figure 1) in this report. Because both 6 and $ dihedral angles can be defined in this compound, it is considered a simple model of a dipeptide even though it does not contain two distinct glycine residues. The usefulness of the potential-annealing method was also examined, using the same model system. All calculations were performed with AMBER version 4.014 slightly modified by us and its all atom force field.17J8 The initial 6 and $values of 1were both set at -180O. This compound was placed in a rectangular box filled with approximately 800 TIP3P water molecules.I6 The system was then minimized until the root mean square value of the potential gradient was below 0.03 kcal mol-' A-'. Following the minimization, normal MD and the potential-scalingMD calculations were performed under periodic boundary conditions. The system was maintained at constant pressure (1 atm) and temperature (300 K) during the runs by Berendsen's weak coupling method.I9 A 9-A nonbonded cutoff was chosen and the SHAKE algorithm15was used to constrain all bond lengths. In the potential-scaling MD, only the bond angle and dihedral potential terms in the potential energy function were reduced by the scaling factor p = 2.0, as mentioned above. The normal MD calculationswere done under the same conditions as in the potential-scaling MD, except that the potential energy was not scaled. At first,the location of the energy minima on the conformational space of 1was explored. Figure 2 shows the +$ potential energy surface calculated with the AMBER potential field. The data points in this figure were generated as follows. Each set of 4 and $ angles of 1was harmonically constrained at 20° intervals with a force constant of 2000 kcal mol-' A-1. All other degrees of freedom were then minimized by using a distance-dependent dielectric constant. The minimizations were performed until the root mean square value of the potential gradient was below 0.03 kcal mol-' A-1 or the step of the minimization r.eached 18 000 cycles. As shown in Figure 2, five minima exist on the conformational energy surface of 1. The minima, A and D, had the same and the lowest energy levels. The energy levels of the
Tsujishita et al.
4418 The Journal of Physical Chemistry, Vol. 97, No. 17, 1993 Potential-Scaled MD
Normal MD
PotentYAnnealing
(a)
PotentinlAnnealing
300 K to 0 K
300 K to 0 K
180 1204
'
'
':
120
60
- r e o - ! . . ...$:' . -180 -120 -60
0
60
120 180
-180 -120 -60
I
0
- . . . 60
120
-
I
i
180
0
3
(a) 20 PS
-60
-120
30
.180 -120 -60
0
60
120 180
-180 -120
40
0
60
120
180
(W 50 ps
(b)
40 Time (ps)
PotentialAnnealing
30
40 Time (ps)
-Potential. Annealing
300KtoOK
180
180
50
300 K to 0 K
120
-E
60
t o z 8 (C)
100 ps
Figure 3. Conformers of 1 appearing during potential-scaled M D and normal M D runs over (a) 20 ps, (b) 50 ps, and (c) 100 ps. and 1c. dihedral angle values of each conformer are plotted. The horizontal and vertical axes in each drawing represent &J and $ dihedral angles, respectively.
-60
-120
-180
other three minima B, C, and E relative to those of A and D were 1.97, 1.92, and 3.43 kJ mol-', respectively. The areas of the minima A and D were much wider than those of the others, and the minimum E had the smallest area. We examined the effectiveness of the potential-scaled MD by observing how efficiently this method could search these minima. Although this potential surface of 1 was evaluated in the system without the explicit solvent, the real (4, $) energy surface in solution is considered to be approximately the same as this one (see Discussion section). Next the potential-scaled MD and the normal MD calculations were performed by using the model system 1 at 300 K. In the potential-scaled MD calculation, the bond angle and dihedral interaction terms in the potential function were scaled down by the factor p = 2. The normal MD calculation was performed under the same condition as in the potential-scaled MD, except for the potential scaling. The results obtained from both calculations are shown in Figure 3. In this figure, the 4 and $ values of the conformers of 1 that appeared in both simulations during the 20-, 50-, and 100-ps runs were plotted. This figure shows how extensively the methods could search the potential energy surface. As indicated by the figure, the potential-scaled MD had already searched a wider area of the potential surface (Figure 2) at 50 ps than the normal M D did, and at 100 ps, the method successfully searched for all of the five energy minima on thesurface.. In contrast, the normal M D could find the minima A, B, and C but failed in searching for the minima D and E in the same period. In addition, most of the conformers of 1 appearing during the potential-scaled M D existed in the low-
-180 90
100
Time (ps)
110
90
100
110
Time (ps)
Figure 4. Dihedral angle transition during the potential-annealing calculations starting from (a) 30 ps and (b) 90 ps structures of the potential-scaled M D simulation.
energy regions of the potential surface and few conformers appeared in the high-energy regions. This indicates the reality of the conformational behavior of the molecule was maintained in the potential-scaled MD. The potential-annealing method, the combination of the potential-scaled M D and simulated annealing, was also tested. Starting with the coordinates and velocities at a certain time in the potential-scaled M D simulation described above, the potential scaling factor p was gradually decreased from 2.0 to 1.0 over 10 ps. This procedure may correspond to changing the temperature of a part of the molecule. At p = 1.0, the potential energy of the selected interactions reverted to the original level. Next the reference temperature was gradually reduced from 300 to 0 K over 10ps. In this test, two potential-annealing calculations were performed, starting with the structures a t 30 and 90 ps in the potential-scaled MD simulation. The system at 30 ps in the potential-scaled MD, whose 4 and $ values were - 6 0 . 8 O and -51S0, respectively, was located in the minimum C on the potential surface. This minimum is a "local" one, which has somewhat higher energy than the lowest energy minima. Figure 4a shows the result of the potential annealing calculated from this 30-ps structure. The 4 value of 1 remained between - 9 O O and-60° through the run. On theother hand, the$valuechanged
Potential-Scaled Molecular Dynamics from -50’ to -180’ at the initial period and the system moved to the edge of the global minimum region A. A successive conformational change occurred near 42 ps, and the system fell deeper into region in A. The final values of 4 and $were - 8 2 . 2 O and 110.lo, respectively, and at last the system reached the deep point of the global minimum A. The other potential-annealing calculation was performed starting with the 90-ps structure in the potential-scaled MD. The 4 and IC, values of this starting structure were -149.4’ and -27.5’, respectively, and the system at this time was located halfway up the energy hill on the potential surface. The result of the calculation is described in Figure 4b. During the potential-annealing process (90-100 ps), each of the 4 and $ dihedral angles rotated in a 360° arc and reverted to its initial value. Then the 4 angle changed to about 80’ a t the end of the annealing process, and finally the system fell to the lowest potential minimum D (4 = 68.3’, $ = -79.5’).
Discussion The purpose of the potential-scaled MD method is to search the conformational spaceefficiently in simulations with the explicit solvent water, maintaining the reality of the simulations with water. The results show that our method could search for all the minima in the conformational space of 1within a shorter period than the normal MD. Only the bond angle and dihedral terms in the potential energy function were reduced to avoid bad influences on the water, and this was sufficient for the method to show its ability in conformational analyses. Moreover, most of the conformers of 1 appearing during the potential-scaled MD existed in the low-energy regions of the potential surface and few conformers appeared in the high-energy regions. This suggests that the topography of the potential surface (i.e., the presence and location of the wells and hills of the energy) was maintained, although the potential energy function was partially modified and the energy level was scaled down in the potential-scaled MD. In spite of the effectiveness of the potential-scaled MD in conformational analyses, the exploration of the global minimum becomes rather difficult if the solute molecule is large and many local minima exist on the conformational space. This is the reason we designed the potential-annealing method. The potentialannealing method can find the low-energy regions in the conformational space “automatically”, as well as the simulated annealing. The utility of the potential annealing was tested with the two conformers appearing in the potential-scaled MD calculation as the initial coordinates. One of these calculations, which started with the conformer located in the “local” energy minimum C, gave the conformer at the lowest potential well A. The other calculation was performed from the conformer halfway up to the potential hill near the minimum C and found the lowest energy minimum D. Both calculations were successful in searching for the lowest energy regions on the potential energy surface. There exist two lowest energy regions, A and D, on the potential surface of 1 (Figure 2). One of the potential-annealing calculations successfully found the minimum A, and the other found the minimum D. The system used here has hundreds of explicit water molecules and includes a huge number of degrees of freedom. Therefore, when an ordinary energy minimizer was applied to these initial conformers, the system could not find the same energy minima (unpublished data). The results suggest the possibilities of the potential-annealing method as a useful tool for searching for low-energy regions in conformational space. Our examinations and discussions up to now have been based on the (4, $) potential energy surface (Figure 2) generated by the energy minimizations of 1. These minimizations were carried out in a pseudovacuum condition, using a distance-dependent dielectric constant to compensate for the lack of explicit water. As mentioned above, however, the presence of explicit solvent can affect the conformational behavior of solute and therefore it is possible that the surface indicated in Figure 2 may differ
The Journal of Physical Chemistry, Vol. 97, No. 17, 1993 4419 from the one in solution. For more accurate discussions, the potential (or free) energy surface of 1 in solution should be calculated. This can be done, for example, by the potential of mean force and free energy perturbation method as Pearlman et al. did,20 but this method requires too much computer cost (a few nanoseconds of M D simulations are needed). In the test system used here, there are no internal hydrogen bonds or strong electrostatic interactions that are greatly affected by the presence of water, since the compound 1 is a short peptide and does not include ionic residues such as glutamic acid or lysine. Therefore it is considered that solvent water causes little effect on the shape of the (4, $) energy surface in this system and that the shape of the (4, $) energy surface in solution is similar to that of the surface in pseudovacuum (Figure 2). The conformers appearing in the 100-ps normal M D simulation (Figure 3) existed a t lowenergy regions of the surface in pseudovacuum and not at highenergy regions. If the topography of the (4, fi) energy surface in solution much differed from that in pseudovacuum, the distribution of the conformers could not be observed as in Figure 3. However, the above description is generally not true in the case of a biomolecule such as a protein or a nucleotide. Such a molecule usually has many internal hydrogen bonds and electrostatic interactions that contribute to the molecular stability. As we mentioned above, these interactions are greatly affected by polar water solvent and therefore the conformational behavior of the molecule will be quite different in solution and in uucuo. It is very difficult to know how the energy hypersurface of the molecule will change when the solvent water exists. Nevertheless, the methods described here are based on molecular dynamics, which is the well-established technique to simulate the real behavior of a system, and therefore they expected to have the ability of searching for low-energy regions on the hypersurface, even though the accurate shape of the surface is unknown. This study proved that the potential-scaled MD method is more effective in conformational analyses than the normal MD. However, a rather long simulation time (about 100 ps) was required, although the model compound used here was small. With a larger system, the method might require more simulation time to display its ability sufficiently. For application to a larger molecular system, further improvements are needed for this method to be more efficient. In the potential-annealing method, the conformational space should be searched as widely as possible in the first (potential-reduced) stage in order to show its usefulness adequately. One approach for the enhancement of the power of scanning the conformational space is the use of a higher value for the potential scaling factor p. From eq 2, a high value of p in the potential-scaled M D corresponds to a high temperature in the high-temperature molecular dynamics. We reported here only the case of p = 2, but it is expected that the potential-scaled M D with a larger p value would search the conformational space more efficiently. The potential-scaled MD with various values is now under study. Another approach considers the solvent effect on conformational changes in solute molecules. Recent work determined the rate constant for the conformational change of a small molecule using Langevin dynamics (LD) simulations with no explicit water.*’ In LD, collisions with the solvent are approximated by the stochastic term in Langevin’s equation of motion, and their magnitude is described by thecollision frequency y. In that work, the relationship between the isomerization rate and the collision frequency y was examined. The results showed that there was an optimal value of y for the maximum transition rate. This was the result of a balance between two solvent effects opposing to each other. At higher y than the optimum, the motion of the solute is rapidly damped by the solvent collisions. At lower 7,energy activation by the solvent is not sufficient to surmount the energy barriers. Our potential-scaled MD and potentialannealing methods include the explicit solvent water molecules,
4420 The Journal of Physical Chemistry, Vol. 97, No. 17, 1993 and this corresponds to the LD with a high y value. The solvent molecules therefore may reduce the effect of the potential-scaled MD by their collisions with the solute. In order to enhance the effectiveness of the potential-scaled MD and potential annealing, some modification of the factors that govern the behavior of the solvent molecules (nonbonded interactions, the weak coupling to the heat bath, etc.) may be needed.
Conclusion An effective conformational search method, potential-scaled MD, was examined in the MD simulationswith the explicit solvent water. This method was tested in a simple model and found to be capable of searching the conformational space accurately and more efficiently than the normal MD. Using the idea of the potential-scaled MD and simulated annealing, we designed the potential-annealing method and examined its usefulness. This method was able to find the lowest energy regions in the conformational space of the model system, starting with the different conformers. These methods need further modifications to enhance their efficiency for a larger biomolecular system, but the results presented here suggest the potency of these methods as useful tools for conformational analyses. Because MD simulations in solution generally yield more accurate results than simulations in vacuo, the potential-scaled MD and potentialannealing methods may be effective for predicting the threedimensional structures of biomolecules. Although this work was performed with a small model compound, the methods proposed here can be expected to show their abilities in simulations of larger systems. The improvement of these methods is under study, and we are also applying.them to a larger peptide system. The results will be reported elsewhere.
Tsujishita et al.
References and Notes (1) Bruccoleri, R. E.; Karplus, M. Biopolymers 1990, 29, 1847-1862. (2) Di Nola, A.; Berendsen, H. J. C.; Edholm, 0.Macromolecules 1984, 17, 20442050. (3) Kirkpatrick, S.; Gelatt, J. C. D.; Vecchi, M. P. Science 1983, 220, 671680. (4) Briinger, A. T.; Kuriyan, J.; Karplus, M. Science 1987, 235, 458460. (5) Briinger, A. T.; Karplus, M.; Petsko, G. A. Acta Crystallogr. 1989, A45, 50-61. (6) Kuriyan, J.; Briinger, A. T.; Karplus, M.; Hendrickson, W. A. Acta Crystallogr. 1989, A45, 396409. (7) Holak, T. A.;Gondol, D.;Otlewski, J.; Wilusz, T. J. Mol. Biol. 1989, 210, 635448.
( 8 ) Nilges, M.; Gronenborn, A. M.; Briinger, A. T.; Clore, G. M. Protein Eng. 1988, 2, 27-38. (9) Zarbock, J.; Oschkinat, H.; Hannappel, E.; Kalbacher, H.; Voelter, W.; Holak, T. A. Biochemistry 1990, 29, 7814-7821. (10) van Gunsteren, W. F.; Berendsen, J. C. Angew. Chem., Int. Ed. Engl. 1990, 29, 992-1023. (11) van Gunsteren, W. F.; Mark, A. E. Eur. J. Eiochem. 1992, 204, 947-96 1. (12) Levitt, M.;Sharon, R. Proc. Natl. Acad.Sci. U.S.A.1988,85,75577561. (13) Mark, A. E.; van Gunstern, W. F.;Berendsen, H. J. C. J. Chem. Phys. 1991, 94, 3808-3816. (14) Pearlman, D. A.; Case, D. A.; Caldwell, J. C.; Seibel, G. L.; Singh, U. C.; Weiner, P.; Kollman, P. A. AMBER 4.0, University of California, San Francisco, 1991. (15) van Gusteren, W. F.; Berendesen, H. J. C. Mol. Phys. 1977, 34, 1311-1327. (16) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (17) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.;Weiner, P. J. Am. Chem. SOC.1984,106,765-784. (18) Weiner, S. J.; Kollman, P. A. J . Comput. Chem. 1986, 7, 230-252. (19) Berendsen, H. J. C.; Postma, J. P. M.;vanGunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984,81, 3684, (20) Pearlman, D. A.; Kollman,P. A. J.Am. Chem. SOC.1991,113,71677177. (21) Loncharich, R.J.; Brooks, B. R.; Pastor, R. W. Eiopolymers 1992, 32, 523-535.
