J. Phys. Chem. B 1998, 102, 6419-6424
6419
Molecular Dynamics Simulations in Aqueous Solution: Application to Free Energy Calculation of Oligopeptides Koichi Tazaki* and Kentaro Shimizu Department of Biotechnology, The UniVersity of Tokyo, 1-1-1 Yayoi, Bunkyo-ku, Tokyo 113-8657, Japan ReceiVed: December 3, 1997; In Final Form: June 2, 1998
We developed a new method to calculate thermodynamic parameters with the Poisson-Boltzmann equation (PBEQ) and molecular dynamics (MD) calculations. This method treats the solvent effect on the solute conformation as a potential of mean force (PMF), so no explicit water is required in the calculation. The usage of the PBEQ enables the reduction of conformational space to be explored, and we used it to calculate the conformational free energy of alanine dipeptide and hydrogen bond forming-deforming (helix-coil) transition free energy of oligopeptides Ace-Ala3-NMe in aqueous solution. We found that (i) the conformational free energy of alanine dipeptide agrees with previous works, (ii) the hydrogen bond formingdeforming transition free energy Ghelix - Gcoil was 1.27 ( 0.27 kcal/mol, and this deviation was small in comparison with the explicit water calculation, (iii) the solvent effect stabilized the hydrogen bond forming conformation, and (iv) electrostatic dipole moment and accessible surface area of polar atoms have an important role for the hydrogen bond forming-deforming transition free energy.
I. Introduction One of the most important areas of computational chemistry is the calculation of thermodynamic parameters, such as free energy surface and specific heat, and the factors determining the stability of R-helices and β-sheets have been of major interest for many years. One of the approaches to understanding these factors is to calculate the thermodynamic parameters of helixcoil transition.1-10 In principle, the thermodynamic parameters are derived from the partition function Z, which is defined as
Z)
∫exp[-E(r1, r2, ..., rN)/kBT]dr1 dr2 ... drN
(1)
where N is the atom number, kB is the Boltzmann constant, T is the temperature, ri is the position of the ith atom, and E is the energy function. In general, the conformational space is too wide for eq 1 to be calculated straightforwardly, so Monte Carlo (MC) simulation and molecular dynamics (MD) simulation, etc., are used to explore the conformational space for the integration of eq 1. The free energy perturbation method (FEP)11,12 is one of the most typical methods that use the MD for conformational exploration. In the case of proteins and nucleic acids, the solvent term in eq 1 must not be omitted because the solvent plays important roles in their structures and functions. But to include the solvent term requires an enormous conformational space so that the naive MC and MD simulations are not valid (Figure 1). Various methods for the exploration have been developed,13-15 but it is difficult to improve the accuracy of the exploration. Improving the accuracy requires either that the size of the conformational space be reduced or that a large amount of statistical data be used for each conformation. The latter is not practical because the conformational space increases exponentially with increasing degree of freedom.
Figure 1. Sampled conformational space by MD, denoted by the hatched area. The area is too small to cover the whole conformational space.
One of the ways to reduce the size of the conformational space is to use the potential of mean force (PMF). The mean force Emean is given as
Emean(T, {r}n) ) -kBT log
∫exp[-E({r}n, ..., {r}N)/kBT]drn+1 ... drN
(2)
where {r}n denotes r1, r2, ..., rn. If {r}n are the coordinates of the protein atoms and rn+1, rn+2, ..., rN are the coordinates of the solvent atoms, this equation represents the integration of the solvent term in eq 1. Therefore, the explicit calculation of the size of solvent-solution interactions is not needed, and the size of the conformational space is considerably reduced. Beglov et al.16 estimated the PMF as an analytical form by using the RISM-HNC integral equations17 and continuum solvent model. In this method, the water molecules that surround the solute were explicitly treated so the conformational space is too wide to explore. The Poisson-Boltzmann equation (PBEQ) was introduced to approximate the PMF of the water. Because the solution of the PBEQ is one of the approximations
S1089-5647(98)00378-2 CCC: $15.00 © 1998 American Chemical Society Published on Web 07/29/1998
6420 J. Phys. Chem. B, Vol. 102, No. 33, 1998 of the PMF, the incorporation of the PBEQ into MD and MC simulations are required. MD requires a force calculation; i.e., the differential of the energy. Niedermeier et al.18 wrote that the differential is the product of the atomic charge and the electrostatic field in aqueous solution. Gilson et al.19,20 reported another method that merges the PBEQ and MD approaches. In these two papers the forces due to the solvent were calculated. But in the methods used in the papers, summation of the forces was not equal to zero because Niedermeier neglected the dielectric boundary force and Gilson et al. took wider grid spacing than the spacing reported in ref 19. In general, MD carried out the temperature scaling, and the nonzero net force makes the scaling difficult. Sharp reported the calculation of the differential using the atomic charge and induced charge (DIEMOND).21 DIEMOND is a useful method, but the induced charge is a conformational dependent term, so the treatment of this charge is a difficult problem. In the paper, the charges were varied through the conformational exploration, but it is hard to know when the charges should be recalculated. The present paper describes a new method that incorporates the solution of the PBEQ into conformational exploration and applies it for the free energy calculation. This method is composed of two phases. One is the MD phase, which uses the DIEMOND method, but we did not change the induced charge. The other is the data analysis phase. In the MD phase, we imposed a harmonic restrain potential into the conformational energy and carried out MD calculations. Because the conformation is restricted, we fixed the induced charges through the simulation. To widely explore the conformational space, we carried out the calculation more than once by altering the harmonic energy. The induced charges were calculated for each simulation. In the data analysis phase, the results of the MD were processed by the weighted histogram analysis method (WHAM).3,7,22-24 We calculated the conformational free energy of alanine dipeptide at first. The alanine dipeptide is small enough to calculate the PMF, and many studies calculated the PMF and free energy of its conformation.15,16,20,21,25-29 So, the calculation of the dipeptide was carried out for the validation of this method. Next, we calculated the hydrogen bond forming-deforming transition free energy of the simplest model helix, Ace-Ala3-NMe.1-4 The model helix contains one hydrogen bond in helix-like conformation, so the hydrogen bond forming-deforming transition indicates the helix-coil transition. We found that (i) the conformational free energy of alanine dipeptide agrees with previous works, (ii) the hydrogen bond forming-deforming transition free energy Ghelix - Gcoil was 1.27 ( 0.27 kcal/mol, and this deviation was small in comparison with the explicit water calculation, (iii) the solvent effect stabilized the hydrogen bond forming conformation, and (iv) electrostatic dipole moment and accessible surface area of polar atoms have an important role for the hydrogen bond forming-deforming transition free energy. II. Method II.1. Molecular Dynamics Simulation with Solution of the Poisson-Boltzmann Equation. In a previous paper30 we developed a method that incorporates the solution of the PBEQ into the Monte Carlo (MC) simulation. That method is summarized as 1. solving the PBEQ numerically by the finite differential method (FDPB) and 2. calculating the “effective charge”, which mimics the spatial distribution of the solution of PBEQ. The sum of the Coulombic energy and the energy due to the solvent, E solv, was given by
Tazaki and Shimizu
E solv )
1 2
∑ i,j
(q ej qi + q ei qj)
1 exp[-κsrij] 4πs
rij
(3)
where q ej and q j are the jth effective and atomic charges, s and ks are the dielectric constant and Debye’s shielding parameter of the solvent, and rij is the distance between the ith and jth atoms. Because the solvent effect was included in the “effective charge” q e and the exponential term, once the PBEQ was solved no more FDPB calculations were required. In reality the “effective charge” is dependent on the molecular conformation, and the “effective charge” was recalculated every several thousand MC steps. The computational time was nevertheless much less than in simulation with FDPB calculations in every MC step. In the previous paper, this method was used to calculate the conformational preferences of the alanine dipeptide and gave results in agreement with the experimental results. But MC calculation is inefficient for large molecules, so application to MD is required. The present paper therefore describes our development of an MD algorithm incorporating the solution of the PBEQ. In the previous method, q e is a conformation-dependent term, so the force due to the solvent and Coulombic interactions, F solv, is given as
)F solv i
1 2
(q ej qi + q ei qj)∇i φ(rij) ∑ i,j 1
(qi∇i q ej + qi∇i q ei )φ(rij) ∑ 2 i,j
(4)
where φ(r) is exp (-ksr)/4πsr. The deviation of the “effective charge” is small, so we assumed the differential of q e to be zero and derived the force on the ith atom, F solv i , as
F solv ) i
1
∑(q ej qi + q ei qj) 2 i*j
1 + κsrij exp[-κsrij] 4πs
r2ij
(ri - rj) (5)
To keep the conformational change small, we added the harmonic conformational restrain energy to the conformational energy, and we did not change the charge in the molecular dynamics simulation. II.2. Conformational Integration. The existence of large energy barriers between two energy minima makes exploration of the conformational space difficult and as a result makes the estimation of the partition function inaccurate. To overcome this problem, we introduce the conformational restrain energy into the conformational energy. Altering the conformational restrain energy, we carried out more than one MD simulation and merged the MD results by using the weighted histogram analysis method (WHAM) developed by Kumar et al.22 The addition of the restrain energy results in the diminution of the energy maxima, so each estimation of the partition function will be more accurate. The altering of the conformational restrain energy will result in the covering of the whole conformational space (Figure 2). In the WHAM conformational density, Ω(ξ, E), where ξ is a reaction coordinate and E is the conformational energy, is given by
Ω(ξ, E) )
∑i Ni(ξ, E)/gi nm exp[-E/kBTm - Em(ξ)/kBTm]/gmZm ∑ m
(6)
MD Simulations in Aqueous Solutions
J. Phys. Chem. B, Vol. 102, No. 33, 1998 6421
Figure 3. Definition of dihedral angles, (φ, ψ), of the alanine dipeptide. Figure 2. Sampled conformational space by each conformation constrained MD, denoted by each ellipse. The restrain energies of each MD are different from each other, so the sampled area of each MD is different from each other. Each area is too small to cover the whole conformational space, but the union of the areas is large enough to almost cover the whole conformational space.
where
Zi )
Ω(ξ, E)exp[-E/kBTi - Ei(ξ)/kBTi] ∑ ξ,E
(7)
and Ti is temperature at the ith simulation, Ni(ξ, E) is the number of observations at (ξ, E) during the ith simulation, ni is the total number of observations, Ei is the conformational restrain energy of the ith simulation, and gi is the correlation time of the ith simulation. We followed Fincham et al.31,32 for the determination of the correlation time. Equations 6 and 7 were iterated to determine Ω(ξ, E) and Zi. The partition function in eq 1 is rearranged by using the conformational density as follows
Z)
Ω(ξ, E) exp[-βE] ∑ ξ,E
(8)
where β is 1/kBT. Various thermodynamic parameters are given in the same way. II.3. Materials and Simulation Procedures. II.3.1. Alanine Dipeptide. To calculate the PMF of the dihedral angles, (φ, ψ), of the alanine dipeptide, we divided the (φ, ψ) plane into 144 windows, i.e., each window covered about a 30° × 30° region. To constrain the conformation near each window, we added the constrain energy, 10 cos(φ - φi) + 10 cos(ψ ψi) kcal/mol, where (φ, ψ) indicated the center position of the each window, to the conformational energy for each simulation. The SHAKE algorithm was used to fix the length of all bonds that included hydrogen. The time step of the simulations was 1 fs. In the first 65 ps, the molecule was gradually heated to 300 K, in the successive 25 ps the temperature was restrained to 300 K, and in the successive 10 ps no temperature scaling procedure was adopted. The conformational sampling was carried out for successive 200 ps, so the total sampling time was 28.8 ns. In this calculation in the solvent, we set the solvent dielectric to 80 and dielectric of the alanine dipeptide to 1.0. II.3.2. Ace-Ala3-NMe. We applied the method to the calculation of the hydrogen bond forming-deforming transition free energy of the Ace-Ala3-NMe depicted in Figure 4. The molecules contain one hydrogen bond in helix-like conformation, so hydrogen bond forming-deforming transition indicates the helix-coil transition. We used as the restrain function in eqs 6 and 7 the harmonic function 1.2(r - ri)2 kcal/mol, where r is the distance between the oxygen in the acetyl group and the hydrogen in the methylamide group. The equilibrium length ri ranged from 1.8 to 13.8 Å at intervals of 1 Å. That is, we simulated the oligopeptides 13 times.
Figure 4. Distance between the oxygen in the acetyl group and the hydrogen in the methylamide group, used as a reaction coordinate.
In the FDPB procedure, the parameters of the solvent were calculated from the pure SPC/E water33 simulation.34 That is, the dielectric constants of the solvent and the peptide were set to 64.1 and 1, and the radius of the water was set to 1.78 Å. We set 0.65 Å for grid spacing. Comparison with the “explicit water model” is undergoing. In the MD procedure we used a force field with the program package with AMBER4.1. “Effective charge” for each simulation was calculated from the initial conformation of each simulation, and the value was fixed through the simulation. The SHAKE algorithm was used to fix the length of all bonds that included hydrogen. The time step of the simulations was 1 fs. In the first 65 ps, the molecule was gradually heated to 300 K, in the successive 25 ps the temperature was restrained to 300 K, and in the successive 10 ps no temperature scaling procedure was adopted. Then the reaction coordinate, the conformation energy, and the conformational constrain energy were recorded every 50 fs for 1 ns in NVE ensemble. That is, there is 100 ps equilibration steps and 1 ns sampling steps in one MD. III. Results and Discussion III.1. Validation by Alanine Dipeptide. To verify our method, we calculated the PMF of the dihedral angles, (φ, ψ), of the alanine dipeptide (Figure 3). The alanine dipeptide is small enough to calculate the PMF, and many studies calculated the PMF and free energy of its conformation.15,16,20,21,25-29 So the calculation of the PMF of the alanine dipeptide is good measure for the justification of our method. It is consensus that the C7eq, (φ, ψ) ) (-70°, 60°), and C7ax, (φ, ψ) ) (70°, -60°) conformations are dominant in the vacuum because of the electrostatic repulsion of the carbonyl oxygens. On the other hand, it is consensus that solvent produces many local minimum in the PMF of the (φ, ψ), and especially stabilizes the RR, (φ, ψ) ) (-60°, -50°), and RL, (φ, ψ) ) (60°, 50°) conformation, which represent the helix-like conformation because of the reduction of repulsion caused by the high dielectric. Previous studies25-29 estimated that solvent reduces the free energy difference between C7 and RR to 1.3 kcal/mol. Figure 5 shows the PMF of the (φ, ψ) for vacuum and aqueous solution. Clearly C7 conformer is dominant in the vacuum (Figure 5a). In contrast to the result of the vacuum, RR and RL conformers are stabilized in aqueous solution (Figure 5b). The free energy differences between various conformations
6422 J. Phys. Chem. B, Vol. 102, No. 33, 1998
Tazaki and Shimizu
Figure 6. “Effective charge” of the oxygen in the acetyl group of each conformational restrain energy. The differences between windows 1 and 13 were large, but the differences between the values for two adjacent windows were small.
TABLE 2: Helix Coil Transition Energy Derived from 10 Different Sets of Initial Conformations
Figure 5. PMF of (φ, ψ) of the alanine dipeptide in a vacuum (a) and in aqueous solution (b). See text for details.
TABLE 1: Free Energy Difference between C7eq, C7ax, rR, rL, and C5a environment
C7eq f RR
C7ax f RL
RR f C5
vacuum aqueous solution Anderson et al.36 Pettitt et al.17 Tobias et al.29
2.9 1.4 1.4 1.3
2.7 0.6 1.1 1.1 0.5
0.0 0.5 0.1 0.2
a All values are shown in kcal/mol. The free energies in solution are similar to those of previous studies.
are shown in Table 1. The energy differences are close to those of the previous studies. These results indicate that our method was a good approximation of the FEP of the solvent effect on the solute conformation. III.2. Hydrogen Bond Forming Free Energy: Application to the Ace-Ala3-NMe. We assumed that the “effective charge” in eq 3 is constant and derived eq 5. Because the “effective charge” depends on the molecular conformation, the differential of the charge is not zero and eq 5 is not exact. However, the error will be small because the molecular conformation is restrained through each MD calculation. The “effective charge” of the oxygen in the acetyl group in each window is shown in Figure 6. The maximum deviation was 20%, but the deviation for two adjacent simulations was always within 5%. This small deviation of the “effective charge” in each MD calculation is a result of the restriction of the molecular conformation. The “effective charge” used in the each window depends on the initial molecular conformation. To determine whether the assumption used in deriving eq 5 from eq 3 was good, we carried out conformational exploration 10 times with different sets of initial conformations. Then the trajectories were analyzed by WHAM, and the helix-coil
run
∆G (kcal/mol)
δ∆G (kcal/mol)
run
∆G (kcal/mol)
δ∆G (kcal/mol)
1 2 3 4 5
1.23 1.54 1.39 1.21 1.19
0.19 0.30 0.18 0.16 0.30
6 7 8 9 10
1.00 1.40 1.09 1.32 1.22
0.29 0.21 0.23 0.22 0.21
transition free energy ∆cfhG was obtained. The helix state was defined as one for which the reaction coordinate is less than 3 Å. The helix-coil transition free energy ∆cfhG was given by
PT (r < 3) ∆cfhG ) Ghelix - Gcoil ) -kBT ln 1 - PT (r < 3)
(9)
where PT(r < 3) denotes the probability that the reaction coordinate is less than 3 Å. In the WHAM, the relative statistical error of the conformational density, δΩ(ξ, E)/Ω(ξ, E), is given by
δΩ(ξ, E) Ω(ξ, E)
-1/2 ∑i g-1 i Ni(ξ, E)]
)[
(10)
and the statistical error of the probability is given by
δ2PT (r < 3) )
δ2Ω(r, E) exp[-2βE]/Z2 ∑ E,r 8 Å while the effect is complicated in the region where r < 8 Å. The effect related to the ASA of the polar atoms increases, but the effect related to the electrostatic dipole moment decreases. These features were consistent with Figure 9. So, we consider that these two factors have an important role on the solvent effect on the conformational preferences.
(1) Yang, A. S.; Honig, B. J. Mol. Biol. 1995, 252, 351. (2) Wang, L.; O’Connell, T.; Tropsha, A.; Hermans, J. Proc. Natl. Acad. Sci. U.S.A. 1995, 92, 10924. (3) Boczko, E. M.; Brooks, C. L., III. J. Phys. Chem. 1993, 97, 4509. (4) Tobias, D. J.; Brooks, C. L., III. Biochemistry 1991, 30, 6059. (5) Hermans, J.; Anderson, A. G.; Yun, R. H. Biochemistry 1992, 31, 5646. (6) Vorobjev, Y. N.; Scheraga, H. A.; Honig, B. J. Phys. Chem. 1995, 99, 7180. (7) Kumar, S.; Payne, P. W.; Va´squez, M. J. Comput. Chem. 1996, 17, 1269. (8) Young, W. S.; Brooks, C. L., III. J. Mol. Biol. 1996, 259, 560. (9) Wang, J.; Purisima, E. O. J. Am. Chem. Soc. 1996, 118, 995. (10) Okamoto, Y.; Hansmann, U. H. E. J. Phys. Chem. 1995, 99, 11276. (11) Singh, U. C.; Brown, F. K.; Bash, P. A.; Kollman, P. A. J. Am. Chem. Soc. 1987, 109, 1607. (12) Bash, P. A.; Singh, U. C.; Lnagridge, R.; Kollman, P. A. Science 1987, 236, 564. (13) Ferguson, D. M.; Pearlman, D. A.; Swope, W. C.; Kollman, P. A. J. Comput. Chem. 1992, 13, 362. (14) Saito, M.; Nakamura, H. J. Comput. Chem. 1990, 11, 76. (15) Hermans, J.; Yun, R. H.; Anderson, A. G. J. Comput. Chem. 1992, 13, 429. (16) Beglov, D.; Roux, B. J. Chem. Phys. 1994, 100, 9050. (17) Pettitt, B. M.; Karplus, M. Chem. Phys. Lett. 1985, 121, 194. (18) Niedermeier, C.; Schulten, K. Mol. Simul. 1992, 8, 361. (19) Gilson, M. K.; Davis, M. E.; Luty, B. A.; McCammon, J. A. J. Phys. Chem. 1993, 97, 3591. (20) Gilson, M. K.; McCammon, J.A.; Madura, J. D. J. Comput. Chem. 1995, 16, 1081. (21) Sharp, K. J. Comput. Chem. 1991, 12, 454. (22) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992, 13, 1011. (23) Kumar, S.; Rosenberg, J. M.; Bouszida, D.; Swendsen, R. H.; Kollman, P. A. J. Comput. Chem. 1995, 16, 1339. (24) Boczko, E. M.; Brooks, C. L., III. Science 1995, 269, 393. (25) Schmidt, A. B.; Fine, R. M. Mol. Sim. 1994, 13, 347. (26) Smith, P. E.; Pettitt, B. M.; Karplusm M. J. Phys. Chem. 1993, 97, 6907. (27) Mezei, M.; Mehrotra, P. L.; Beveridge, D. L. J. Am. Chem. Soc. 1985, 107, 2239. (28) Brady, J.; Karplus, M. J. Am. Chem. Soc. 1985, 107, 6103. (29) Tobias, D. J.; Brooks, C. L., III. J. Phys. Chem. 1992, 96, 3864. (30) Tazaki, K.; Doi, J. J. Phys. Chem, 1996, 100, 14520. (31) Fincham, D.; Quirke, N.; Tildersley, D. J. J. Chem. Phys. 1986, 84, 4535. (32) Allen, M. P.; Tildersley, D. J. Computer Simulation of Liquids, 2nd ed.; Oxford University Press: New York, 1990; Chapter 6. (33) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (34) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (35) Zimm, B. H.; Bragg, J. K. J. Chem. Phys. 1959, 31, 526. (36) Anderson, A. G.; Hermans, J. Proteins 1988, 3, 262.
