Replica Exchange Monte Carlo Simulations for Folding of Di-Block Polyampholyte
Yuri YAMADA, Yosuke UEDA and Yosuke KATAOKA
Complex systems have a huge number of local minima in the phase space. Especially at low temperature, the molecular configuration tends to be trapped in the local minima and the comprehensive investigation of the phase space is obstructed. Therefore, the search for the lowest energy configurations of such systems is difficult using the conventional Monte Carlo (MC) and/or molecular dynamics (MD) techniques. So far, various devices have been suggested to overcome this problem: the cluster algorithm [1, 2], the multicanonical method [3 - 6], the simulated tempering [7 - 9], and so on. Replica Exchange (RE) algorithm [10 - 13], that can be regarded as a parallelized version of the simulated tempering, has been also widely used and been confirmed as an effective method in the search for the global minimum of such complex systems. Further information for such algorithms for MC and/or MD is found in Ref. .
Polyampholyte (PA) [5, 6, 15 - 18], a chain-like molecular model, is a typical case of the complex systems. We simulated the folding of the PA using the Replica Exchange Monte Carlo (REMC) method.
2 Replica Exchange Method
The RE method assumes a compound system consisting of plural replicas. The replicas are equivalent and independent of each other. While the replicas have common initial configuration, their temperatures are different from each other. The replicas are simulated simultaneously and independently. The exchanges are performed during certain periods of the simulation steps, and the exchanging pair of replicas is chosen randomly among all the replicas. By the exchanges, the configuration trapped on the local minimum will surmount the potential barrier, and will arrive at the global minimum.
To keep the equilibrium of the overall system through the RE process, it is necessary that the system satisfies the detailed balance condition. For the exchanging configuration of the m-th and n-th replicas, the detailed balance condition is given with the transition matrix W :
Here, the symbols Xm, bm, Em, and Tm are the configuration, the inverse temperature, the energy, and the temperature of m-th replica, respectively, and kB is the Boltzmann constant. Thus, the exchange probability is expressed as:
Another virtue of the RE is easiness for parallel computing. The replicas are simulated independently except for the period of the exchange. Accordingly, the RE calculation can be fastened ideally in parallel with the number of replicas.
3 Model: Di-block Polyampholyte
Polyampholyte (PA) is a chain-like molecule model that consists of charged monomers, like peptides and proteins [5, 6, 15 - 18]. We treat "di-block PA" which has an elementary distribution of the charges; the halves of monomers are posi-tively charged on one side of the chain and the other halves are negatively charged on the other side . This is illustrated in Figure 1. We assume the intramolecular interaction of the di-block PA, U, as the sum of the Coulombic interaction Uc between each pair of monomers, the elastic interaction Ue between the adjoining monomers, and the softcore interaction Us between each pair of monomers excluding the adjoining pairs:
Here, qi is the electric charge of i-th monomer, e0 is the vacuum permittivity, rij is the distance between i-th and j-th monomers, k and n are the elastic parameters, re is the equilibrium distance between adjoining monomers, e and s are softcore parameters. We use the reduced units for the length rr, the energy Er, and the temperature Tr as:
where e is the elementary electric charge.
Figure 1. The illustration of the di-block polyampholyte.
We performed the REMC calculations for folding of the di-block PA. The conventional MC (Simulated Annealing: SA or SAMC) calculations were also performed for a comparison. The calculation procedures are stated below, and the parameters and conditions are shown in Table 1.
The MC program was coded by Fortran 77, and the calculations performed with Alpha-Linux on non-parallelized system: because some variant calculations were performed simultaneously. We used the multiplicative method for random number generation.
Table 1. The parameters and conditions of the REMC and the SAMC calculations.
a We consider N Metropolis trial moves in one MC step.b The RE trials are performed periodically every these MC steps. The number of MC steps is increased in the latter periods (on lower temperatures) by 5.0 × 105 steps.c The products of the RE trials and the MC steps on the corresponding periods of the REMC.
|number of monomers: N||number of replicas||number of RE trials||MC stepsa / RE trialb||MC steps for the SAc|
|60||19||100||1.0 × 106 ~ 3.0 × 106||1.0 × 108 ~ 3.0 × 108|
|q / e||re / r0||e / E0||n||k / E0r0-4|
|±1||1||1||4||1.6 × 104|
Figure 2. The illustration of our procedure of the REMC.
REMC. The number of monomers of the di-block PA (N) is 60. The system contains one PA molecule in a vacuum. We consider the system to be a replica, and we assume 19 replicas to be an overall system. The initial configurations of all replicas are the same: a straight chain. The upper limit of the temperature is 100 E0/kB, and we set the temperatures in 5 E0/kB steps into the 19 replicas. Concretely, the temperatures of the replicas are 100, 95, 90... 15, and 10 E0/kB. For such compound system, the conventional Metropolis MC calculations are performed independently on each replica, and the RE trials are done periodically. We perform 100 RE trials in one period of calculation. After the termination of the first period, we adopt the configuration of a replica with the lowest temperature (T = 10 E0/kB) as the initial configuration of the next period, and the temperatures of all replicas are multiplied by 1/10. Thus, the second period calculation is performed on the temperatures of 10, 9.5, 9.0... 1.5, and 1.0 E0/kB. We repeat such procedure to the fifth period. Consequently, the lower limit of the temperature is 0.001 E0/kB. Such procedure is illustrated in Figure 2.
SAMC. Conditions of SAMC calculations are nearly the same as those of the REMC. Yet, of course, the replicas are not assumed and the temperature of the system is annealed sequentially: as 100, 95, 90... 15, 10, 9.5, 9.0... 0.0015, and 0.001 E0/kB.
5 Results and Discussions
Figure 3 indicates the final snapshots of the calculations. We obtained the lowest energy configuration of our model to be the stretched double-helical structure. On the other hand, the sphere-like helix configurations, trapped in the local minima, were also obtained. We call the former Stretched helix (Figure 3(a)), and the latter Spherical helix (Figure 3(b)).
In Figures 4, 5 we show the temperature dependences of the averaged intramolecular interaction, <U>, and the averaged distance between the center of mass of the PA molecule and each monomer, <r>, respectively. In the low temperatures, the Stretched helix obtains a lower <U> than the Spherical helix, and their difference appears at the temperature ca. Tr <= 0.3. The <r> curves also divide around this temperature, moreover the curve of the Stretched helix has a minimum at Tr ~ 0.2. It proves that the PA molecule surmounts the potential barrier around the local minima and tends toward the global minimum (Stretched helix).
We performed 6 calculations for both REMC and SAMC respectively. The Stretched helix appeared not only in the RE but also in the SA, and the RE was hardly better than the SA for obtaining the probability of the Stretched helix: 3/6 in the RE and 2/6 in the SA. We conceive this result to be caused by the shape of the potential surface; if the potential barrier surrounding the local minimum is low, the local-minimum state easily surmounts the barrier without regard to application of the RE. It is considered that the potential function and its parameters determine the height of the potential barrier.
Shimizu, Hiwatari, and coworkers have performed multicanonical MD calculations for the folding transition of the di-block PA (N = 20, 30, 40, and 60) [5, 6], and concluded that the Stretched helix is the lowest energy configuration. Our model PA was somewhat different from theirs; the distances between adjoining monomers were kept constant in their model. Anyway, it is significant that there is agreement on the lowest energy configuration of the both models.
Figure 3. The final snapshots of the calculations at Tr = 0.001. (a): The lowest energy configuration. (b): An example of the local-minima configuration. These two configurations were obtained by both RE and SA calculations.
Figure 4. The plots of the averaged intramolecular interaction <U> as a function of the temperature Tr. This plot is the magnification of the low-temperature region.
Figure 5. The plots of the averaged distance between the center of mass of the PA and each monomer <r> as a function of the temperature Tr.
We calculated the folding of di-block PA using the REMC and SAMC simulations. It was significant that the lowest energy configuration was predicted as the double-helical structure. As a result of some RE and SA calculations, the double-helical structure appeared not only in the RE but also in the SA. The probability of obtaining the double helix by using the RE was not quite as good as that with the SA. It is considered that this result is caused by the potential function and its parameters.
[ 1] R. H. Swendsen and J-S. Wang, Phys. Rev. Lett., 58, 86 (1987).
[ 2] U. Wolff, Phys. Rev. Lett., 62, 361 (1989).
[ 3] B. A. Berg and T. Neuhaus, Phys. Lett. B, 267, 249 (1991).
[ 4] B. A. Berg and T. Neuhaus, Phys Rev. Lett, 68, 9 (1992).
[ 5] H. Shimizu, K. Uehara, K. Yamamoto, and Y. Hiwatari, Mol. Sim., 22, 285 (1999).
[ 6] A. Baumketner, H. Shimizu, M. Isobe, and Y. Hiwatari, J. Phys.: Condens. Matter, 13, 10279 (2001).
[ 7] A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov, and P. N. Vorontsov-Velyaminov, J. Chem. Phys, 96, 1776 (1992).
[ 8] E. Marinari and G. Parisi, Europhys. Lett., 19, 451 (1992).
[ 9] W. Kerler and P. Rehberg, Phys. Rev. E, 50, 4220 (1994).
 K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn., 65, 1604 (1996).
 K. Hukushima, H. Takayama, and K. Nemoto, Int. J. Mod. Phys. C, 7, 337 (1996).
 M. C. Tesi, E. J. J. van Rensburg, E. Orlandini, and S. G. Whittington, J. Stat. Phys., 82, 155 (1996).
 D. A.Kofke, J. Chem. Phys., 117, 6911 (2002).
 A. Mitsutake, Y. Sugita, and Y. Okamoto, Biopolymers (Pept. Sci.), 60, 96 (2001).
 P. G. Higgs and J-F. Joanny, J. Chem. Phys., 94, 1543 (1991).
 Y. Kantor, H. Li, and M. Kardar, Phys. Rev. Lett., 69, 61 (1992).
 Y. Kantor, M. Kardar, and H. Li, Phys Rev. E, 49, 1383 (1994).
 M. Tanaka, A. Y. Grosberg, V. S. Pande, and T. Tanaka, Phys. Rev. E, 56, 5798 (1997).