Multicanonical molecular dynamics simulations combined with Metadynamics for the free energy landscape of a biomolecular system with high energy barriers

Chemical Physics Letters 501 (2011) 598–602

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Multicanonical molecular dynamics simulations combined with Metadynamics for the free energy landscape of a biomolecular system with high energy barriers Yasushige Yonezawa ⇑, Hiromitsu Shimoyama, Haruki Nakamura Institute for Protein Research, Osaka University, 3-2 Yamadaoka, Suita, Osaka 565-0871, Japan

a r t i c l e

i n f o

Article history: Received 28 October 2010 In final form 17 November 2010 Available online 19 November 2010

a b s t r a c t Multicanonical molecular-dynamics (McMD) simulation and Metadynamics (MetaD) are useful for obtaining the free-energies, and can be mutually complementary. We combined McMD with MetaD, and applied it to the conformational free energy calculations of a proline dipeptide. First, MetaD was performed along the dihedral angle at the prolyl bond and we obtained a coarse biasing potential. After adding the biasing potential to the dihedral angle potential energy, we conducted McMD with the modified potential energy. Enhanced sampling was achieved for all degrees-of-freedom, and the sampling of the dihedral angle space was facilitated. After reweighting, we obtained an accurate free energy landscape. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Molecular dynamics (MD) and Monte Carlo simulations are indispensable methods that are widely used for studies of condensed matter, such as bio-molecules. Rugged free energy landscapes are observed in biological molecules, in particular, where the stable phase spaces are separated by high free energy barriers as in the case of the cis–trans isomerization of proline residues. Therefore, to elucidate those physiochemical aspects, it is necessary to determine the free energy landscapes accurately. Conventional canonical molecular simulations frequently fail to yield accurate free energies, because of the ruggedness and limited computational time. Consequently, to date, several algorithms have been proposed for obtaining the free energies. Although most of the algorithms are based on MD or Monte Carlo simulations, an alternative method, the reaction path Hamiltonian superposition approach proposed by Strodel and Wales, is based on geometry optimization that provides very fast sampling of transition sate regions [1]. The multicanonical algorithm was first proposed by Berg and Neuhaus in Monte Carlo simulations [2]. It was later extended to molecular dynamic simulations (McMD) [3,4], and greatly improves the sampling efficiency of the simulations. McMD simulations are a generalized Boltzmann ensemble algorithm that uses the temperature as a general reaction coordinate. Once the multicanonical weight is determined accurately, the reweighting procedures produce the canonical distribution at an arbitrary temperature. McMD simulations have been effectively applied to several biological targets [5,6]. Despite the wide energy searching

⇑ Corresponding author. Fax: +81 66879 8636. E-mail address: [email protected] (Y. Yonezawa). 0009-2614/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2010.11.061

capability, McMD simulations do not provide satisfactory sampling of the phase spaces around the transition states with free energies that are much higher than the stable states. Although Ono and coworkers resolved this problem by combining the McMD simulations with the umbrella sampling method, in a technique referred to as multicanonical WHAM [7], their method, which performs many McMD simulations in each umbrella sampling window, is difficult and very time consuming. Metadynamics (MetaD), an algorithm developed recently by Laio and Parrinello, effectively enhances the sampling and evaluation of the free energy landscapes [8]. In MetaD, the system is simulated using conventional MD simulations, and a historydependent repulsive biasing potential is added to the phase space of a few collective variables. The biasing potential is usually a Gaussian form with fixed parameters, and it discourages the system from approaching previously visited configurational space. These parameters, which are system dependent, should be chosen ad hoc and carefully tuned to improve the performance. The history-dependent biasing potential rapidly fills the free energy minima and enables the efficient exploration of the phase space. Once the free energy surfaces are flattened, it is possible to walk in the phase space randomly. Therefore, MetaD has been combined with several other methods such as umbrella sampling [9,10], parallel tempering [11], and the Wang–Landau method [12]. Although the combinations are useful, a reliable estimation method for the free energy errors is still needed [9,13]. In this report, we introduce a method that combines McMD with MetaD. We then use it to investigate the free energy landscape of a bio-molecule with high energy barriers, such as transition states. We show that the two methods, McMD and MetaD, can be complementary to each other, and achieve effective sampling along a reaction coordinate across high free energy barriers. Although an improved MetaD algorithm has been proposed to

599

Y. Yonezawa et al. / Chemical Physics Letters 501 (2011) 598–602

address the convergence issue [13], herein we utilized the rapid evaluation feature for the free energies in the original MetaD algorithm. We chose a proline dipeptide as a biological model system. The proline dipeptide has two stable states, trans and cis, where the trans state is more stable than the cis states by about 1.5 kcal/mol in water. Since the two states are separated by free energy barriers larger than 14 kcal/mol, it is impossible to overcome them using a conventional MD simulation at room temperature. We applied MetaD to a reaction coordinate with respect to the cis–trans isomerization and obtained a biasing potential, which reproduces the approximate free energy along the reaction coordinate. In this way, the application of the biasing potential to the reaction coordinate is expected to reduce the free energy barriers and to facilitate the following McMD sampling along the reaction coordinate. We show that the McMD simulation with the MetaD modification provides much more effective sampling than that without the modification. Consequently, we demonstrated that the combination of MetaD and McMD is an effective method for obtaining accurate free energies of bio-molecules with high energy barriers.

Me1

We modeled the simulation system by immersing N-acetyl-N’methylprolineamide (Ace-Pro-NMe) in the center of a pre-equilibrated TIP3P water [14] cubic box with 23.7 Å edges, and removing the overlapping water molecules within 3.0 Å from the dipeptide. The number of water molecules in the system was 446. The center-of-mass of the Ace-Pro-NMe was fixed to avoid drifting during the simulations. We minimized the energy of the system by 300 steps of the conjugate gradient and 100 steps of the steepest descent methods. Next, we gradually equilibrated the system from the low temperature of 3 K to the high temperature of the production runs using NVT molecular dynamics simulations with the Hoover–Evans method [15] over 10 ps. The SHAKE algorithm [16] was applied to all covalent bonds associated with hydrogen atoms. Amber force field (parm96) [17] was used to model the atomic interactions. We used the particle mesh Ewald (PME) method [18] for the long range electrostatic interactions. The a parameter of the PME method was set to 0.35 Å1, and the real space cutoff distance was set to 8 Å for all runs. The mesh size was set to 24  24  24, thereby ensuring the grid density of 1 Å for the system with sufficient accuracy from the Ewald method. In addition, we used the group base cutoff scheme. A cutoff of 8 Å was used for the treatment of the short range vdW interactions. We integrated the equations of motion using the velocity Verlet integration algorithm with a time step of 1 fs. The final conformation of the peptide was the trans state with x0 nearly equal 180°, where x0 is the dihedral angle of O1–C1–N– Cd, as shown in Figure 1. Note that our definition of x0 is different from the standard dihedral angle x of peptide backbones. Fischer et al. used the angle f for prolyl cis–trans isomerization, instead of x, to decouple the change of the isomerization from that of the pyramidalization [19]. Similarly, x0 was used to perform the uncoupling and to create a reaction coordinate for the umbrella sampling simulations [20]. 2.2. Metadynamics molecular dynamics simulation According to the original MetaD report [8], we used the standard biasing potential management procedure in the MetaD simulations. The biasing potential has the following form,

H N

H

C2

ψ

Cα

Cδ

O2

H

H

Cγ H

Cβ

H

H

Figure 1. Definition of the dihedral angle and atoms for the calculations of the proline dipeptide, Ace-Pro-NMe.

U bias ðx ; tÞ ¼

2.1. Preparation of the proline dipeptide system

H2

ω'

0

2. Method

Me2

C1

O1

X ti

ðx0  x0 ðt i ÞÞ2 h exp  2w2

! ð1Þ

Each Gaussian is centered at scheduled updating points, x0 (ti), along the dihedral angle x0 at time ti. Here, h is the height and w is the width. The updating period, the height, and the width are important parameters used to determine the accuracy and effectiveness of MetaD. We set h and w to 0.25 kcal/mol and 0.5 rad1, respectively. The biasing potential was added every 1000 steps. Effective and rapid convergence was achieved at 400 K in the MetaD simulations. The production run was performed at 400 K. We did not pay particular attention when determining the MetaD parameters because the following McMD simulations could correct the incomplete biasing potential. Consequently, the modified potential energy of the system is described as:

Emod ¼ E þ U bias ðx0 Þ

ð2Þ

Therein, Emod denotes the modified potential energy, consisting of the original potential energy, E, and the additional energy along the dihedral angle obtained by using MetaD. In addition, Ubias(x0 ) stands for the converged biasing potential, Ubias(x0 ,t), in Eq. (1). 2.3. Multicanonical molecular dynamics simulations We used the force-biased Multicanonical molecular dynamics (FBMcMD) simulation method [21,22] to evaluate the free energy landscape. The FBMcMD algorithm was described previously [21– 23], and so we will briefly summarize it. We generated the FBMcMD ensemble by performing constant-temperature MD simulations at an arbitrarily chosen multicanonical temperature, Tmc = 1/kBbmc, with force scaling, as presented below:

dpi ¼ mðEÞfi dt @ a ðEÞ mðEÞ ¼ mc @E 1 amc ðEÞ ¼ E þ ln PðEÞ bmc

ð3:1Þ ð3:2Þ ð3:3Þ

Therein, bmc represents the inverse multicanonical temperature, kB signifies the Boltzmann’s constant, and pi and fi stand for the momentum and the force of the i-th atom, respectively. In addition, amc(T) signifies the multicanonical energy and m(E) represents the force scaling factor. Since m(E) has not been given a priori, it should be determined using the following iterative scheme:

mkþ1 ðEÞ ¼ mk ðEÞ þ

1 @ ln Pk ðEÞ bmc @E

ð4Þ

600

Y. Yonezawa et al. / Chemical Physics Letters 501 (2011) 598–602

Here, Pk(E) is the probability distribution of potential energy from the k-th iterative run of the FBMcMD. Furthermore, m(E) is related to the density of states, X(E), through the following equation:

1

XðEÞ

¼ ebmc amc ðEÞ

ð5Þ

Once m(T) converges, the system can achieve a random walk on the potential energy space. We set the temperature range for the FBMcMD simulation to 300–700 K, so that the peptide could sample various structures at these temperatures. We performed two FBMcMD simulations: one simulation used the original potential energy, E; the other one used the modified potential energy, Emod. We set the multicanonical temperature to 200 K. We then obtained flat energy distributions covering the temperature range from 300–700 K, as shown in Figure 2, executed production runs for 2  107 steps, and stored snapshots every 1 ps. All simulations were performed using the myPresto molecular dynamics computing program [24]. 2.4. Reweighting in the combined method The potential of mean force at a desired temperature, with the original potential energy from the combined method, is reproduced using reweighting techniques, as described in the following. We first perform the standard multicanonical reweighting technique to a desired temperature, T0, using Eq. (6):

ebmc amc ðEmod Þ eb0 Emod PðEmod ; bo Þ ¼ Z0

ð6:1Þ

where

Z0 ¼

Z

dEmod ebmc amc ðEmod Þ eb0 Emod

ð6:2Þ

In these equations, P(Emod, b0) represents the probability distribution produced from the modified potential energy Emod at T0 = 1/ kBb0. T0 stands for the reweighting temperature, and Z0 is the partition function. We subtract the biasing potential, Ubias(x0 ), from the potential of mean force, F(Emod, b0), calculated from Eq. (7). Here, we can ignore the contribution of Z0, because it does not depend on energy and temperature.

FðE; b0 Þ ¼ FðEmod ; b0 Þ  U bias ðx0 Þ

ð7:1Þ

where

FðEmod ; b0 Þ ¼ kT 0 ln PðEmod ; b0 Þ

ð7:2Þ

3. Results and discussion 3.1. Biasing potential along the reaction coordinate From the MetaD simulations, we obtained the biasing potential along the reaction coordinates, as depicted in Figure 3. The minus sign of the biasing potential is similar to the potential of mean force (PMF) along the dihedral angle. Consequently, adding the biasing potential to the dihedral angle potential energy can reduce the barrier height. To verify the reduction, we added the biasing potential to the PMF obtained by the umbrella sampling MD simulations along the dihedral angle conjugated with the umbrella integration method [25]. The reduced potential is presented in Figure 3. The fact that the maximum value is 3.8 kcal/mol indicates the sufficient reduction of the large PMF barriers. The reduction is expected to facilitate the sampling ability of the next McMD. The incomplete reduction is not a problem, because the dihedral angle space with the modified potential is expected to be well sampled and corrected by the following McMD simulations, which is the remarkable advantage of the current method.

3.2. Modified system calculated using Multicanonical molecular dynamics simulation Using the system with the modified dihedral angle potential, we performed McMD simulations with 2  106 steps. The time course of the dihedral angle x0 from the simulations is shown in Figure 4. Moreover, we performed McMD simulations with the temperature range of 300–1000 K to investigate higher temperature effects. The simulations started from 40 different initial conformations of the trans state. As a result, 39 of the McMD simulations maintained the trans state during the simulations, while one McMD simulation only once transformed from the trans to cis states. The time course of the dihedral angle x0 of the transformed simulation, as well as the time course of the McMD simulation with temperatures of 300–700 K, is presented in Figure 4. The figures show that the McMD simulation with the modified potential was not trapped in the vicinity of the initial state, and x0 moved around the whole area of the dihedral angle space, although almost all of the conventional McMD simulations with the original potential were trapped in the trans state. Furthermore, even when the conventional McMD simulation covered a much higher temperature, 1000 K, the energy barriers between the trans and cis states were rarely overcome. Thus, McMD with the modified po-

Finally, we obtain the potential of mean force, F(E, b0), as derived from the original potential energy E at the temperature T0.

Figure 2. The logarithm of the probability as a function of the total potential energy obtained by using the McMD simulations. The range of energies corresponds to the temperatures from 300–700 K.

Figure 3. Solid line: the biasing potential obtained from the MetaD simulations. Dashed line: the reduced potential produced by adding the biasing potential to the potential of mean force obtained from the umbrella sampling simulations. The cis and trans states correspond to the vicinities of 0° and ±180° in the dihedral angle, respectively.

Y. Yonezawa et al. / Chemical Physics Letters 501 (2011) 598–602

601

Figure 4. Time courses of the dihedral angle x0 from the McMD simulation with the bias potential over the temperature range of 300–700 K (upper panel); from the McMD simulation over the temperature range of 300–1000 K, transforming the trans to cis states (middle panel); and from the McMD simulation with the modified potential (lower panel).

tential from MetaD facilitated the structure sampling quit effectively. The potential of mean force (PMF) of the dipeptide, as functions of x0 and w reweighted at 300 K by the McMD simulations with a temperature range of 300–1000 K, is presented in Figure 5. The same PMF by the McMD simulations with the modified potential is presented in Figure 6. The one-dimensional projection of the PMF, as a function of x0 , is also displayed in Figure 6. Figure 5 shows that McMD without the modified potential searches only the restricted narrow area, because of the large barriers and the poor statistics. In contrast, as shown in the upper panel of Figure 6, McMD simulations with the modified potential can search a much wider area than that without the potential modification, revealing that the biasing potential facilitates the sampling ability of the McMD simulations remarkably. Moreover, the

Figure 5. The conformational potential of mean forces of the prolyl dipeptide from the 300–1000 K McMD simulation reweighting to 300 K. The simulation transforms the trans to cis states. Unit is in kcal/mol.

Figure 6. Upper panel: (a) the conformational potential of mean force of the prolyl dipeptide from the McMD simulations combined with MetaD reweighted to 300 K, in which the biasing potential is subtracted. Unit is in kcal/mol. Lower panel: (b) one-dimensional projection of the potential of mean force, as a function of x0 .

lower panel of Figure 6 shows that the trans state is more stable, by 1.2 kcal/mol, than the cis state, and that the free energy barriers separating the two states are about 15 kcal/mol. Both are in good accordance with previously reported values [26]. Although the MetaD provides the coarse biasing potentials with respect to the true free energy barriers, the results reveal that the McMD simulations corrected the error due to the incomplete flatness of the energy barriers by the coarse biasing potentials, and they reproduced the accurate free energy landscape. It must be emphasized that no particular care is necessary to determine the MetaD parameters in the current method, while much more efforts is usually required in the original MetaD method. Ono and co-workers combined the umbrella sampling method and McMD simulations to achieve effective sampling in the cis– trans isomerization of a proline peptide [7]. Their method and the current method are logically equivalent. However, many McMD simulations, corresponding to the number of the umbrella sampling windows, are necessary with their method. In contrast, the current method requires only a single McMD simulation. Namely, the biasing potential in the current method is regarded as a global umbrella potential along the reaction coordinate. A similar method was reported by Kannan and Zacharias, in which a peptide backbone biasing potential combined with replica exchange simulation provided the enhanced peptide backbone transitions [27]. The biasing potential produced by the umbrella sampling method requires several umbrella simulations, which are necessary to build a continuous biasing potential along the reaction coordinate and numerous replicas for the replica ex-

602

Y. Yonezawa et al. / Chemical Physics Letters 501 (2011) 598–602

change MD simulations. In contrast, our current method requires only two sequential simulations. Babin and co-workers reported the free energy landscape of a small peptide using MetaD, and they performed the correction using additional umbrella sampling simulations [9]. This combination method resembles the current method. However, the current method not only corrects the phase space sampled by MetaD, but also achieves enhanced sampling of all of the degrees of freedom, whereas their method corrects a limited number of degrees of freedom. It has been pointed out that Amber force field (parm96) does not provide a very realistic free energy barrier of the cis–trans isomerization for proline [20,28]. Its value is several kcal/mol smaller than that obtained from experiments [28]. For obtaining the accurate free energy barriers strongly relating with their electronic structures, the hybrid QM/MM MD simulation [20] should be required. Consequently, our results show that when appropriate degrees of freedom, which cannot be overcome by McMD simulations, are known, then MetaD can rapidly flatten the degrees of freedom by the biasing potentials. Then, using the potential modified by MetaD, the McMD simulations significantly improve the sampling ability and quality. These results show that the current combined method is useful for the free energy calculations of biological molecules. 4. Conclusion To investigate the accurate free energy of bio-molecules efficiently, we proposed a method combining Multicanonical molecular dynamics (McMD) with Metadynamics (MetaD) simulations. Although McMD offers an important advantage to the sampling capability independent of the conformational reaction coordinates, McMD simulations present a disadvantage in sampling around transition states that are much higher than the stable phase spaces. In contrast, the original MetaD has the advantage of rapid evaluation of the free energy along a reaction coordinate, while the reaction coordinates are usually limited to small numbers and the conversion is generally difficult. Our combined method complementarily incorporates the two advantages. If the reaction coordinates around the transition states are known, then we can apply MetaD to the reaction coordinates to obtain a biasing potential. Adding the biasing potential to the reaction coordinate potential reduces the free energy barriers. Subsequently, McMD simulations can achieve effective sampling around the transition states including other degrees of freedom, because of the barrier reduction.

Using the combined method, we conducted the conformational free energy evaluations of a proline dipeptide, which undergoes transitions between the cis and trans states and has much higher energies than those at the stable states. By reweighting the combined method to 300 K, we obtained an accurate free energy landscape, which demonstrates that the combined method is useful in general for accurate and rapid studies of the free energy of biomolecules with high energy barriers. Acknowledgement This research was supported by Research and Development of the Next-Generation Integrated Simulation of Living Matter, a part of the Development and Use of the Next-Generation Supercomputer Project of the Ministry of Education, Culture, Sports, Science and Technology (MEXT). References [1] B. Strodel, D.J. Wales, Chem. Phys. Lett. 466 (2008) 105. [2] B. Berg, T. Neuhaus, Phys. Rev. Lett. 68 (1) (1991) 9. [3] H.E.H. Ulrich, Y. Okamoto, F. Eisenmenger, Chem. Phys. Lett. 259 (3–4) (1996) 321. [4] N. Nakajima, H. Nakamura, A. Kidera, J. Phys. Chem. B 101 (1997) 817. [5] U.H.E. Hansmann, Y. Okamoto, J. Comput. Chem. 14 (11) (1993) 1333. [6] N. Nakajima, J. Higo, A. Kidera, H. Nakamura, J. Mol. Biol. 296 (2000) 197. [7] S. Ono, N. Nakajima, J. Higo, H. Nakamura, Chem. Phys. Lett. 312 (1999) 247. [8] A. Laio, M. Parrinell, Proc. Natl. Acad. Sci. USA 99 (20) (2002) 12562. [9] V. Babin, C. Roland, T.A. Darden, C. Sagui, J. Chem. Phys. 125 (2006) 204909. [10] B. Ensing, A. Laio, M. Parrinello, M. Klein, J. Phys. Chem. B 109 (2005) 6676. [11] G. Bussi, F.L. Gervasio, A. Laio, M. Parrinello, J. Am. Chem. Soc. 128 (2006) 13435. [12] D. Min, Y. Liu, I. Carbone, W. Yang, J. Chem. Phys. 126 (2007) 194104. [13] A. Barducci, G. Bussi, M. Parrinello, Phys. Rev. Lett. 100 (2008) 020603. [14] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [15] D.J. Evans, W.G. Hoover, B.H. Failor, B. Moran, A.J.C. Ladd, Phys. Rev. A 28 (1983) 1016–1046. [16] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [17] W.D. Cornell et al., J. Am. Chem. Soc. 117 (1995) 5179. [18] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, J. Chem. Phys. 103 (1995) 8577. [19] S. Fischer, R.L. Dunbrack Jr., M. Karplus, J. Am. Chem. Soc. 116 (1994) 11931. [20] Y. Yonezawa, K. Nakata, K. Sakakura, T. Takada, H. Nakamura, J. Am. Chem. Soc. 131 (12) (2009) 4535. [21] J.G. Kim, Y. Fukunishi, A. Kidera, H. Nakamura, Phys. Rev. E 68 (2003) 21110. [22] J.G. Kim, Y. Fukunishi, A. Kidera, H. Nakamura, Phys. Rev. E 70 (2004) 57103. [23] Y. Watanabe, J.G. Kim, Y. Fukunishi, H. Nakamura, Chem. Phys. Lett. 400 (2004) 258. [24] Y. Fukunishi, Y. Mikami, H. Nakamura, J. Phys. Chem. B 107 (2003) 13201. [25] J. Kastner, W. Thiel, J. Chem. Phys. 123 (2005) 144104. [26] C. Melis, G. Bussi, C.R. Lummis, C. Molteni, J. Phys. Chem. B 113 (2009) 12148. [27] S. Kannan, M. Zacharias, Proteins 66 (2007) 697. [28] U. Doshi, D. Hamelberg, J. Phys. Chem. B 113 (2009) 16590.