Replica-Permutation Method with the Suwa-Todo Algorithm beyond the Replica-Exchange Method
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] O c t Replica-Permutation Method with the Suwa-Todo Algorithmbeyond the Replica-Exchange Method
Satoru G. Itoh ∗ and Hisashi Okumura † Department of Theoretical and Computational Molecular ScienceInstitute for Molecular ScienceOkazaki, Aichi 444-8585, JapanDepartment of Structural Molecular ScienceThe Graduate University for Advanced StudiesOkazaki, Aichi 444-8585, Japan
Abstract
We propose a new method for molecular dynamics and Monte Carlo simulations, which is referredto as the replica-permutation method (RPM), to realize more efficient sampling than the replica-exchange method (REM). In RPM not only exchanges between two replicas but also permutationsamong more than two replicas are performed. Furthermore, instead of the Metropolis algorithm,the Suwa-Todo algorithm is employed for replica-permutation trials to minimize its rejection ratio.We applied RPM to particles in a double-well potential energy, Met-enkephalin in vacuum, and aC-peptide analog of ribonuclease A in explicit water. For a comparison purposes, replica-exchangemolecular dynamics simulations were also performed. As a result, RPM sampled not only thetemperature space but also the conformational space more efficiently than REM for all systems.From our simulations of C-peptide, we obtained the α -helix structure with salt-bridges betweenGly2 and Arg10 which is known in experiments. Calculating its free-energy landscape, the foldingpathway was revealed from an extended structure to the α -helix structure with the salt-bridges.We found that the folding pathway consists of the two steps: The first step is the “salt-bridgeformation step”, and the second step is the “ α -helix formation step”. ∗ Electronic address: [email protected] † Electronic address: [email protected] . INTRODUCTION In recent years, generalized-ensemble algorithms are frequently employed to studybiomolecular systems (for reviews, see, e.g., Refs. [1, 2]). This is because it is difficultto obtain sufficient sampling in the conformational space of these systems by conventionalcanonical-ensemble simulations [3–9]. The canonical simulations tend to get trapped in alocal-minimum free-energy state of the biomolecular systems.The replica-exchange method (REM) [10, 11] is one of the most well-known methodsamong the generalized-ensemble algorithms. Non-interacting copies (replicas) of a targetsystem are employed in REM. Different temperatures are assigned to the replicas. By ex-changing the temperatures between the replicas, random walks of the replicas in the temper-ature space are realized. Accordingly, the simulation can escape from local-minimum states.It is easier to perform replica-exchange molecular dynamics (MD) or Monte Carlo (MC) sim-ulations than to perform multicanonical MD [12, 13] or MC simulations [14, 15], althoughthe multicanonical algorithm is also one of the most well-known generalized-ensemble algo-rithms. In the multicanonical and similar algorithms [16–25] non-Boltzmann weight factorsare used as the weight factors. These non-Boltzmann weight factors are not a priori knownand have to be determined by iterative procedures. In REM, however, the usual Boltzmannweight factor is employed for each replica. Therefore, there is no necessity to perform theprocedures of obtaining the weight factor.In REM the Metropolis algorithm [3] is utilized to exchange the temperatures betweenthe replicas. The Metropolis algorithm is a Markov chain Monte Carlo (MCMC) methodand widely employed to obtain a required ensemble. The Metropolis algorithm is designedso as to satisfy the detailed balance condition, which is a sufficient condition to performstate transitions based on MCMC. Trials of such state transitions are accepted or rejectedstochastically, and their rejection ratio increases when there is a large difference betweenthe probability distributions. In order to minimize the rejection ratio, new algorithm wasproposed recently by Suwa and Todo [26].In the Suwa-Todo algorithm the detailed balance condition is not imposed. They intro-duced a graphical procedure called weight allocation instead of solving the detailed balancecondition algebraically. By minimizing the rejection ratio with the weight allocation, this al-gorithm realizes efficient sampling of states. For a system which has only two states for each2article like the Ising model, however, this algorithm is exactly the same as the Metropolisalgorithm [26].To take advantage of the Suwa-Todo algorithm, one might consider to exchangetemperatures in REM by this algorithm. For example, let us consider the follow-ing exchange of the temperatures T and T between the replicas, Replica 1 andReplica 2: (Replica 1 at T ; Replica 2 at T ) → (Replica 1 at T ; Replica 2 at T ). In thiscase, the number of all combinations of the replicas and temperatures is only two,(Replica 1 at T ; Replica 2 at T ) and (Replica 1 at T ; Replica 2 at T ). For this exchange,therefore, the Suwa-Todo algorithm is exactly the same as the Metropolis algorithm. Ingeneral, because temperatures are exchanged between two replicas as in this example, REMis not able to take advantage of the Suwa-Todo algorithm.In order to solve this difficulty, we propose a new generalize-ensemble algorithm, whichwe refer to as the replica-permutation method (RPM). In this method temperatures arepermuted among more than two replicas by the Suwa-Todo algorithm. We can utilize theSuwa-Todo algorithm because the number of all combinations of the replicas and tempera-tures is larger than two.We first apply this new method to particles in a double-well potential energy. For acomparison purposes, RPM with the Metropolis algorithm, and REM are also performed. Asthe second application of RPM, we employ Met-enkephalin in vacuum. This penta-peptideis often used as a test system to see the usefulness of new algorithms [27–30]. Furthermore,RPM is applied to a C-peptide analog of ribonuclease A in explicit water, which is knownto form a α -helix structure [31–35], to see its sampling efficiency for a larger biomolecularsystem. The results of the second and third applications are compared with those of REM.It is considered that the α -helix structure of the C-peptide analog is stabilized by saltbridges (SBs) between Gly2 and Arg10 [34]. We discuss the role of the SBs for the α -helixstructure. Furthermore, a folding pathway based on a free-energy landscape is presentedfrom our simulations.In Section II we describe formulation of RPM. The Suwa-Todo algorithm and the graph-ical procedure, weight allocation, for MCMC are also introduced in this section. We presentthe details of our simulations in Section III. The results are shown in Section IV. Section Vis devoted to conclusions. 3 I. METHODSA. Markov Chain Monte Carlo Method with the Weight Allocation
We first describe usual MCMC method with the Metropolis and Suwa-Todo algorithms.Let us consider that a system has n states and that the weight of state i is expressed by w i ( i = 1 , · · · , n ). In MCMC the weight w i is given by w i = n X j =1 w j P ( j → i ) , (1)where P ( j → i ) is the transition probability from state j to state i . By defining the amountof stochastic flow from state j to state i as v ( j → i ) ≡ w j P ( j → i ) , (2)Eq. (1) can be rewritten as w i = n X j =1 v ( j → i ) . (3)This amount of stochastic flow v ( j → i ) also satisfies w j = n X i =1 v ( j → i ) , (4)because n X i =1 P ( j → i ) = 1 . (5)From Eqs. (3) and (4), the global balance equation is derived: n X i =1 v ( i → j ) = n X i =1 v ( j → i ) . (6)By performing state transitions with v ( j → i ) satisfying this equation, a required ensembleis obtained.In the Metropolis algorithm the detailed balance condition, which is a sufficient conditionfor Eq. (6), is employed to obtain the required ensemble. Here, the detailed balance conditionis described in terms of the amount of stochastic flow as follows: v ( i → j ) = v ( j → i ) . (7)4 ( j → i ) is given by the following equation so as to satisfy Eq. (7): v ( j → i ) = 1 n − w j , w i ] , j = i , (8)where the coefficient 1 / ( n −
1) comes from a random selection of state i from ( n − j .When v ( j → i ) satisfies Eqs. (3) and (4), Eq. (6) is automatically fulfilled . Therefore,we focus on Eqs. (3) and (4) after this. These two equations can be understood visuallyby the weight allocation as follows (see Fig. 1(a)): Equation (3) is regarded as filling a boxwhose size is w i by blocks whose sizes are v ( j → i ) ( j = 1 , · · · , n ) without any vacantspace. Equation (4) is regarded as dividing a block whose size is w j into smaller blockswhose sizes are v ( j → i ) ( i = 1 , · · · , n ). In order to satisfy both equations simultaneouslyfor any i and j , therefore, we prepare blocks w j and boxes w i ( i, j = 1 , · · · , n ). By usingblocks v ( j → i ) divided from w j , all w i are filled without any space. Figure 1(a) shows thisweight allocation of Metropolis algorithm for a system which has four states ( n = 4). Theblock size of v ( j → i ) is calculated from Eq. (8). Red frame blocks represent rejected flow v ( i → i ) ( i = 1 , · · · , P i v ( i → i ) / P i w i .A new algorithm is proposed recently by Suwa and Todo through the weight allocationto minimize the average rejection ratio for state transitions in MCMC. We refer to thisalgorithm as the Suwa-Todo algorithm. In this algorithm v ( j → i ) satisfies Eqs. (3) and(4) without imposing the detailed balance condition in Eq. (7). In the Suwa-Todo weightallocation, w j is divided and w i is filled as follows (see Fig. 1(b)):(i) Choose the state which has the maximum weight. If two or more states have themaximum weight, one of them is chosen. Here, we assume that w is the maximumweight without loss of generality.(ii) Box w is filled by block w ( v (1 → w still remains after filling box w ,try to fill the next box w ( v (1 → w become 0 ( v (1 → , · · · , v (1 → k )).(iii) By using block w , fill boxes in turns from the last partially filled box at Step (ii)( v (2 → k ) , · · · , v (2 → l )). This procedure is repeated for the blocks w , · · · , w n ( v (3 → l ) , · · · ). 5iv) Once all the boxes except w are saturated, box w ( · · · , v ( n → n = 4. In this figure the average rejectionratio P i v ( i → i ) / P i w i is 0. In this Suwa-Todo algorithm, v ( j → i ) is given by v ( j → i ) = max [0 , min [∆ ji , w j + w i − ∆ ji , w j , w i ]] , (9)where ∆ ji ≡ S j − S i − + w , (10)and S i ≡ i X j =1 w j , S ≡ S n . (11)By performing state transitions based on the amount of stochastic flow v ( j → i ) or thetransition probability P ( j → i ) = v ( j → i ) /w j , the required ensemble is obtained.Regarding the rejection ratio, from Eq. (9), v ( i → i ) = max [0 , w − S n ] , i = 1 , , i ≥ . (12)Therefore, the rejection ratio becomes 0, if w ≤ S n . (13)It means a reject-free MC simulation is realized if w is less than or equal to the half of thetotal weight S n = P nj =1 w j . This condition can be fulfilled in most cases. B. Replica-Exchange Method with the Metropolis Algorithm
We consider a system of N atoms with their coordinate vectors and momentum vectorsdenoted by q ≡ { q , · · · , q N } and p ≡ { p , · · · , p N } , respectively. The Hamiltonian H instate x ≡ ( q, p ) is given by the sum of the kinetic energy K and potential energy V : H ( x ) = K ( p ) + V ( q ) . (14)In the canonical ensemble at temperature T , each state x is weighted by the Boltzmannfactor: W B ( x ) = e − βH ( x ) , (15)6here β = 1 /k B T ( k B is the Boltzmann constant).Let us suppose that there are M non-interacting copies (or replicas) of the original systemin the canonical ensemble at M different temperatures T m ( m = 1 , · · · , M ). In REM thereplicas are arranged so that there would be always exactly one replica at each temperature.In other words, there is a one-to-one correspondence between the replicas and temperatures.Therefore, the label i ( i = 1 , · · · , M ) for the replicas is a permutation of the label m ( m = 1 , · · · , M ) for the temperatures, and vice versa: i = i ( m ) ≡ f ( m ) ,m = m ( i ) ≡ f − ( i ) , (16)where f ( m ) is a permutation function of m and f − ( i ) is its inverse.Let X α = n x [ i (1)]1 , · · · , x [ i ( M )] M o = n x [1] m (1) , · · · , x [ M ] m ( M ) o stands for a “state” in REM. Here,the superscript i and the subscript m in x [ i ] m are labels of the replicas and temperatures,respectively. All possible combinations between the replicas and temperatures are labeledby the subscript α . The state X α is specified by the M sets of coordinates q [ i ] and momenta p [ i ] of N atoms in replica i at temperature T m : x [ i ] m ≡ (cid:0) q [ i ] , p [ i ] (cid:1) m . (17)Because the replicas are non-interacting, the weight factor w R ( X α ) for the state X α is givenby the product of Boltzmann factors for each replica i (or at temperature T m ): w R ( X α ) = M Y i =1 exp n − β m ( i ) H (cid:16) x [ i ] m ( i ) (cid:17)o = M Y m =1 exp (cid:8) − β m H (cid:0) x [ i ( m )] m (cid:1)(cid:9) , (18)where i ( m ) and m ( i ) are the permutation functions in Eq. (16).We now consider exchanging a pair of replicas in REM. Suppose we exchange replicas j and k which are at temperatures T m and T n , respectively ( j = i ( m ), k = i ( n )): X α = (cid:8) · · · , x [ j ] m , · · · , x [ k ] n , · · · (cid:9) −→ X β = (cid:8) · · · , x [ k ] m , · · · , x [ j ] n , · · · (cid:9) . (19)From Eq. (8) the amount of stochastic flow v ( X α → X β ) for this replica exchange is givenby v ( X α → X β ) = C min [ w R ( X α ) , w R ( X β )] , (20)and the transition probability P ( X α → X β ) is expressed by P ( X α → X β ) = C min (cid:20) , w R ( X β ) w R ( X α ) (cid:21) , (21)7here C = 1 / ( M C ) if replicas j and k were selected randomly among M replicas. Here, M C is the number of 2-combinations from M elements: M C = M ! / { ( M − } . In theusual REM, replica exchanges are tried only between neighboring two temperatures. In thiscase, C = 1 because M = 2. C. Replica-Permutation Method with the Suwa-Todo Algorithm
We consider to perform a replica permutation among all M replicas as a generalizationof a replica exchange in Eq. (19): X α = n x [ i (1)]1 , · · · , x [ i ( M )] M o −→ X β = n x [ j (1)]1 , · · · , x [ j ( M )] M o , (22)where both i and j are permutation functions and i = j . Not only an exchange between tworeplicas, but also a permutation among more than two replicas are allowed in this method.Note that the number of all possible combinations between the replicas and temperatures is M !. Therefore, the index α of X α takes a value between 1 and M !.In the Metropolis algorithm, the transition probability P ( X α → X β ) for this replica per-mutation is also given by Eq. (21). The replicas are allowed to transit to non-neighboringtemperatures, but P ( X α → X β ) takes a quite small value for such replica permutation.Accordingly, most of such replica permutations are rejected. The number of replica per-mutations in which any replica does not transit to non-neighboring temperatures is givenby [ M ] X n =1 M − n C n , (23)where (cid:20) M (cid:21) = M , for even number of M , M − , for odd number of M . (24)On the other hand, the number of all replica-permutation candidates is M ! − (cid:18)P [ M ] n =1 M − n C n (cid:19) / ( M ! − M , because [ M ] X n =1 M − n C n < M X n =0 M C n = 2 M ≪ M ! . (25)8o avoid this rejection problem, we apply the Suwa-Todo algorithm to the replica per-mutations. As in Sec. II A, we assume that w R ( X ) is the maximum weight without loss ofgenerality. The amount of stochastic flow v ( X α → X β ) is determined by the weight alloca-tion in the same way also as in Sec. II A only by replacing the weight w i to w R ( X α ). FromEqs. (9), (10), and (11), v ( X α → X β ) is given by v ( X α → X β ) = max [0 , min [∆ αβ , w R ( X α ) + w R ( X β ) − ∆ αβ , w R ( X α ) , w R ( X β )]] , (26)where ∆ αβ ≡ S α − S β − + w R ( X ) , (27)and S α ≡ α X β =1 w R ( X β ) , S ≡ S M ! . (28)If w R ( X γ ) ( γ = 1) is the maximum weight more generally, Eqs. (27) and (28) are modifiedas follows: ∆ αβ ≡ S α − S β − + w R ( X γ ) , (29)and S α ≡ α X β = γ w R ( X β ) , for α ≥ γ , M ! X β = γ w R ( X β ) + α X β =1 w R ( X β ) , for α < γ , (30) S ≡ S M ! . A replica-permutation simulation with the Suwa-Todo algorithm is performed as follows:Step 1: The label α ( α = 1 , · · · , M !) of X α is assigned to all combinations between the replicasand temperatures. Table I shows an example for M = 3 (three replicas).Step 2: For each replica, a canonical MD or MC simulation at the assigned temperature iscarried out simultaneously and independently for a certain steps.Step 3: A replica-permutation trial is performed as follows: First, each weight is obtainedby Eq. (18), and the maximum weight w R ( X γ ) is determined. Next, we calculatethe amount of stochastic flow v ( X α → X β ) in Eq. (26) and the transition probability P ( X α → X β ) = v ( X α → X β ) /w R ( X α ) for β = 1 , · · · , M !. Finally, transition fromState X α to State X β is accepted stochastically with the probability P ( X α → X β ).9epeating Step 2 and Step 3, we can carry out the replica permutation MD or MC simulation.Figure 2 shows an example of time series of temperatures in RPM. This method realizesnot only minimization of the rejection ratio but also transitions of the replicas to non-neighboring temperatures.The number of combinations between the replicas and temperatures increases in propor-tion to M !. For a large number of replicas, we can divide all replicas and temperatures intosubsets to decrease the number of combinations. Although such a division is not necessary,three to eight replicas are appropriate in each subset. As an example, let us consider thatthe total number of the replicas is eight and that they are divided into two subsets whichhave four replicas. In this case, the replica-permutation simulation is performed as follows:Step 1: Let us suppose that the temperatures are assigned to the replicas at an initial state as Replica 1 at T Replica 2 at T ...Replica 8 at T . They are divided into two subsets:
Replica 1 at T Replica 2 at T Replica 3 at T Replica 4 at T and Replica 5 at T Replica 6 at T Replica 7 at T Replica 8 at T . All combinations between the replicas and temperatures for the former subset andthose for the latter subset are labeled by X α and X α ( α = 1 , · · · , Replica 3 at T Replica 4 at T Replica 5 at T Replica 6 at T and Replica 7 at T Replica 8 at T Replica 1 at T Replica 2 at T . The labels for the former subset and the latter subset are X α and X α , respectively.10tep 2: For each replica, a canonical MD or MC simulation at the assigned temperature iscarried out simultaneously and independently for a certain steps.Step 3: At an odd number of trail time, replica permutations for X α and X α are carried out.That is, four replicas are permutated among four corresponding temperatures in eachsubset. At an even number of trial time, replica permutations for X α and X α areperformed.Step 4: As a result of the replica permutation for X α , let us suppose that the combination ofthe replicas and temperatures is changed as Replica 1 at T Replica 2 at T Replica 3 at T Replica 4 at T → Replica 1 at T Replica 3 at T Replica 2 at T Replica 4 at T . (31)Due to this permutation, the combination in X α is automatically changed without areplica permutation for this subset: Replica 3 at T Replica 4 at T Replica 5 at T Replica 6 at T by Eq. (31) −−−−−−−−−→ Replica 2 at T Replica 4 at T Replica 5 at T Replica 6 at T . (32)To avoid re-labeling for X α , we perform the next permutation regarding Replica 2as Replica 3. Although the combination of the temperatures and replicas may bechanged by the last replica permutation, we do not need to relabel the combinationsas in Table I by replacing the replicas in this way.Repeating Step 2 to Step 4, we continue the replica permutation MD or MC simulation. D. Reweighting Techniques
The results obtained from RPM can be analyzed by the reweighting techniques as inREM [36, 37]. Let us suppose that we have carried out a RPM simulation with M replicasand M different temperatures T m ( m = 1 , · · · , M ).11or appropriate reaction coordinates ξ and ξ , the canonical probability distribution P T ( ξ , ξ ) at any temperature T can be calculated from P T ( ξ , ξ ) = X E M X m =1 ( g m ) − N m ( E ; ξ , ξ ) e − βEM X m =1 ( g m ) − n m e f m − β m E , (33)and e − f m = X ξ ,ξ P T m ( ξ , ξ ) . (34)Here, g m = 1 + 2 τ m , τ m is the integrated autocorrelation time at temperature T m , N m ( E ; ξ , ξ ) is the histogram of the potential energy and the reaction coordinates ξ and ξ at T m , and n m is the total number of samples obtained at T m . Note that this probabil-ity distribution is not normalized. Equations (33) and (34) are solved self-consistently byiteration. For biomolecular systems the quantity g m can safely be set to be a constant inthe reweighting formulas [37], and so we set g m = 1 throughout the analyses in the presentwork. These equations can be easily generalized to any reaction coordinates ( ξ , ξ , · · · ).From the canonical probability distribution P T ( ξ , ξ ) in Eq. (33), the expectation valueof a physical quantity A at any temperature T is given by h A i T = X ξ ,ξ A ( ξ , ξ ) P T ( ξ , ξ ) X ξ ,ξ P T ( ξ , ξ ) . (35)We can also calculate the free energy (or, the potential of mean force) as a function of thereaction coordinates ξ and ξ at any temperature T from F T ( ξ , ξ ) = − k B T ln P T ( ξ , ξ ) . (36) III. COMPUTATIONAL DETAILSA. Asymmetric Double-Well Potential Energy
In order to verify that RPM gives correct ensembles, we first applied this method toa simple system. The system has 100 non-interacting particles with a one-dimensional12symmetric double-well potential energy. The potential energy V ( q ) at a coordinate q isdefined by V ( q ) = (cid:0) ( q + 1) − (cid:1) (cid:0) ( q − − . (cid:1) kcal / (mol · ˚A ) . (37)To see usefulness of the Suwa-Todo algorithm for the replica-permutation method, we per-formed replica-permutation MD simulations with both Metropolis and Suwa-Todo algo-rithms. From now on, we will call the replica-permutation method with the Suwa-Todoalgorithm RPM and that with the Metropolis algorithm M-RPM. The MD versions of the for-mer and latter are called RPMD and M-RPMD, respectively. Conventional replica-exchangeMD (REMD) simulations were also carried out for a comparison purposes. We prepared 40different initial conditions (ICs) for each method. The MD simulation from each IC wasperformed for 10.0 ns including equilibration run for 1.0 ns. The time step was taken to be1.0 fs. The mass of each particle was 1.0 a.u. Six replicas were used, and temperatures wereset at 200 K, 235 K, 275 K, 325 K, 380 K, and 450 K. The temperatures were controlled bythe Gaussian constraint method [4, 5] to avoid the problem of non-ergodicity in the Nos´e-Hoover thermostat [6–8] for the non-interacting particle system. The trajectory data werestored every 10.0 fs. Replica permutations or exchanges were tried every 1.0 ps. B. Met-Enkephalin in Vacuum
To see the usefulness of RPM for biomolecular systems, we next employed a Met-enkephalin molecule in vacuum as a test system. The results were compared with thoseobtained from a REMD simulation. The AMBER parm99SB force field [38, 39] was used.The N-terminus and the C-terminus of Met-enkephalin were blocked by the acetyl groupand the N-methyl group, respectively. Therefore, the amino-acid sequence is Ace-YGGFM-Nme. The SHAKE algorithm [40] was employed to constrain bond lengths with the hydrogenatoms during our simulations. The temperature was controlled by the Nos´e-Hoover thermo-stat [6–8]. The time step was taken to be 1.0 fs. The initial conformation was an extendedstructure.The number of replicas was 12, and temperatures were 200 K, 230 K, 265 K, 300 K,340 K, 385 K, 435 K, 490 K, 555 K, 635 K, 720 K, and 820 K. For PRM, we divided the12 replicas into two subsets, which had six replicas and six temperatures. The RPMD andREMD simulations were performed for 50.0 ns per replica including equilibration run for 1.013s. The trajectory data were stored every 100 fs. Replica permutations or exchanges weretried every 1.0 ps.
C. C-peptide in Explicit Water
In order to demonstrate sampling efficiency of RPM for a lager system, we performeda RPMD simulation of a C-peptide analog in explicit water solvent [32]. This peptide isknown to form an α -helix structure at a lower temperature than 318 K at pH 5.2 [32,35]. Because the charges at the peptide termini affect helix stability [31], we blocked thetermini of C-peptide by neutral Nme and Ace. Therefore, the amino-acid sequence is Ace-AETAAAKFLRAHA-Nme. The histidine residue was protonated to conform our simulationconditions to the experimental pH. A REMD simulation was also carried out. The number ofwater molecules was 1800, and two chlorine ions were added as counter ions. The AMBERparm99SB [38, 39] was used. The model for the water molecules was the TIP3P rigid-body model [41]. Temperature was controlled by the Nos´e-Hoover thermostat [6–8]. TheSHAKE algorithm [40] was employed to constrain bond lengths with the hydrogen atoms ofC-peptide and to fix the water molecule structures during our simulations. This system wasput in a cubic unit cell with the side length of 38.6 ˚A with the periodic boundary conditions.The cutoff distance for the Lennard-Jones potential energy was 12.0 ˚A. The electrostaticpotential energy was calculated by the Ewald method [42]. The multiple-time-step method[43] was employed in our MD simulations. For interactions between the C-peptide atomsand those between the C-peptide atoms and the water atoms, the time step was taken tobe 1.0 fs. The time step of 4.0 fs was used for interactions between the water atoms. Theinitial conformation was an extended structure.The numbers of replicas in the RPMD and REMD simulations were 24, and temperatureswere 281 K, 285 K, 289 K, 294 K, 299 K, 304 K, 309 K, 314 K, 320 K, 326 K, 332 K, 338 K,344 K, 351 K, 358 K, 365 K, 372 K, 380 K, 388 K, 396 K, 405 K, 414 K, 423 K, and 433 K.These replicas were divided into four subsets in RPM. Each subset had six replicas and sixtemperatures. The RPMD and REMD simulations were performed for 40.0 ns per replicaincluding equilibration run for 4.0 ns. Namely, the production run of each simulation wascarried out for 864.0 ns in total. The trajectory data were stored every 400 fs. Trials ofreplica permutations and exchanges were performed every 4.0 ps.14 V. RESULTS AND DISCUSSIONA. Asymmetric Double-Well Potential Energy
We first show that RPM (and M-RPM) gives correct probability distributions for theasymmetric double-well potential energy. The probability distributions P ( q ) in the RPMD,M-RPMD, and REMD simulations are presented in Fig. 3. Here, P ( q ) for each methodwas obtained from the 40 simulations starting from the different ICs. The bin size ∆q for P ( q ) was taken to be 0.05 ˚A. The errors were estimated as the standard deviations of the40 simulations. To see the accuracy of the simulation results, exact probability distributions P exact ( q ) are also illustrated as the solid lines in the figure. P exact ( q ) at a temperature T was calculated numerically by P exact ( q ) = C DW Z q + ∆q q − ∆q dq ′ exp ( − β V ( q ′ )) , (38)where C DW is the normalization constant: C DW = (cid:16)R ∞−∞ dqP exact ( q ) (cid:17) − . The integration inthis equation was computed by the Simpson’s method. As shown in Fig. 3, all P ( q ) showgood agreement with P exact ( q ) at all temperatures. Correct probability distributions areobtained by not only REM but also RPM (and M-RPM).The transition ratios of the replicas from a temperature to another temperature in eachmethod are shown in Fig. 4(a). The transition ratio of the replicas is defined here as aprobability with which a replica at the temperature is transferred to another temperature.RPM has the largest transition ratios at all temperatures among all methods. On the otherhand, those of M-RPM are extremely small (the values are 0.003–0.004). This is becausemost of the replica-permutation trials were rejected. To realize efficient replica permutations,therefore, it is essential to employ the Suwa-Todo algorithm. The total numbers of tunnelingtimes of all replicas during the simulations are listed in Table II. Here, when a replica makes around trip between the lowest and highest temperatures, it is counted as one tunneling. Thenumber of the tunneling times is a useful information to see how efficiently the simulationsamples the temperature space. The value for each method in the table was obtained bytaking an average of the 40 simulations’ results. The number of the tunneling times in RPMwas 1.7 times larger than that in REM. It means that RPM realizes efficient sampling inthe temperature space. On the other hand, M-RPM did not have such a large number of15he tunneling times. This is because its transition ratios were very small.To estimate the convergence speed to the exact probability distributions, we examined thetime series of the average deviations of probability distributions P i,j particle ( q ; t ) from P exact ( q ),where P i,j particle ( q ; t ) is a probability distribution obtained from a single MD simulation ofparticle i from IC j accumulated until time t . The average deviation ∆ q ( t ) at the time t isdefined by ∆ q ( t ) = 1 N bin N IC N particle N bin X k =1 N IC X j =1 N particle X i =1 (cid:12)(cid:12) P i,j particle ( q k ; t ) − P exact ( q k ) (cid:12)(cid:12) , (39)where N bin is the number of the bins, N IC is the number of ICs, N particle is the number ofthe particles in a single simulation, and q k is the coordinate at bin k . 120 bins ( N bin = 120)ranging between q k = − q k = 3 were taken into account. Figure 5 shows ∆ q ( t )calculated from each method. Convergence of M-RPM is the slowest due to the low transitionratios of the replicas. RPM shows slightly faster convergence than REM at the lowesttemperature ( T = 200 K) although convergence at highest temperature ( T = 450 K) isalmost the same. It is important to increase the sampling efficiency at a low temperaturebecause the conformational sampling at a high temperature is originally easier than at a lowtemperature. Therefore, RPM realizes efficient sampling not only in the temperature spacebut also in the coordinate space at a low temperature. It is again because the transitionratios of the replicas by RPM are higher than those by REM. B. Met-Enkephalin in Vacuum
The replica-transition ratios and the total numbers of tunneling times for Met-enkephalinare shown in Fig. 4(b) and Table II, respectively. The transition ratios in RPM are largerthan those of REM at all temperatures. RPM has also larger tunneling times than REM.Thus, efficient sampling in the temperature space was realized by RPM for this biomolecularsystem, too.Its local-minimum free-energy structures in vacuum have been reported for various forcefields although those for the AMBER parm99SB force field [39] have yet to be reported.For example, the global-minimum and local-minimum structures for the CHARMM22 forcefield [44] are shown in Fig. 6. To obtain global-minimum and local-minimum structures inthe AMBER parm99SB force field and to see the sampling efficiency of RPM, we illustrate16ree-energy landscapes at 200 K in Fig. 7. The abscissa and ordinate are the root-mean-square deviation (RMSD) with respect to the structure in Fig 6(a) (RMSD ) and that withrespect to the structure in Fig 6(b) (RMSD ), respectively. Here, the RMSD is defined byRMSD = min s n X j (cid:0) q j − q j (cid:1) , (40)where n is the number of the backbone atoms in Met-enkephalin, { q j } are the coordinatesin the reference conformation, and the minimization is over the rigid translations and rigidrotations for the coordinates of the conformation { q j } with respect to the center of geometry.The free-energy landscapes in Figs. 7(a) and 7(b) were calculated from Eq. (36) with thereweighting techniques. These landscapes show good agreement with each other, and havefive local-minimum free-energy states. The five local-minimum states are labeled as A toE, as shown in the figures. The free-energy landscapes in Figs. 7(c) and 7(d) were obtainedfrom the raw histogram without the reweighting techniques. In the landscape of REM inFig. 7(d), States D and E are not observed although these states are observed in RPM.In order to discuss sampling efficiency at the lowest temperature more quantitatively, wecounted the numbers of visiting times in each state, as listed in Table III. Here, when areplica visited a local-minimum state at the lowest temperature after the replica had visitedanother state at the lowest temperature, it is counted as one visit. The errors were estimatedby the jackknife method [45–47] in which the production run was divided into 20 segments.We regarded the regions presented in Table IV as those for the local-minimum states. Asshown in Table III, REM did not visit State D and State E at the lowest temperature. Thisis the reason why these states were not observed in Fig. 7(d). On the other hand, RPMvisited all states, and the numbers of visiting times in RPM are larger than REM for allstates. RPM thus samples the conformational space more efficiently than REM at the lowesttemperature.The representative conformations at the local-minimum free-energy states are also shownin Fig. 7. The structural features are as follows: (State A) The structure in State A is theglobal-minimum free-energy structure for the AMBER parm99SB force field and similar tothat for the CHARMM22 force field in Fig. 6(a). This structure has two hydrogen bondsbetween NH of Gly2 and CO of Phe4 and between CO of Gly2 and NH of Phe4. Thehydroxy group of the Tyr1 side-chain is close to CO of Gly3. (State B) The structure in17tate B is almost the same as the structure in Fig. 6(b). Two hydrogen bonds are formedbetween NH of Gly2 and CO of Met5 and between CO of Gly2 and NH of Phe4. However,this structure does not have the hydrogen bond between CO of Gly2 and NH of Met5 whichexists in Fig. 6(b). Distance between the hydroxy group of Tyr1 and CO of Gly3 is small,as in State A. (State C) There are two hydrogen bonds between CO of Tyr1 and NH ofPhe4 and between CO of Tyr1 and NH of Met5 in State C. (State D) Two hydrogen bondsbetween CO of Tyr1 and NH of Gly3 and between NH of Gly2 and CO of Met5 are formedin State D. As for the structures in States C and D, the Tyr1 hydroxy group is not close toany backbone CO. (State E) The structure in State E has a hydrogen bond between NH ofGly2 and CO of Phe4. This structure also has a small distance between the Tyr1 hydroxygroup and CO of Met5. C. C-peptide in Explicit Water
The transition ratios of the replicas and the total numbers of tunneling times for C-peptide in explicit water are shown in Fig. 4(c) and Table II, respectively. The transitionratios in RPM are larger than those in REM at all temperatures, again. The tunneling timesof RPM was about 2.1 times larger than that of REM. The time series of the temperaturesof Replica 1, Replica 9, and Replica 17 are shown in Fig. 8. Figure 8 actually shows morefrequent tunneling in RPM than in REM. In REM 6 replicas had never made a round tripbetween the lowest and highest temperatures during the simulation. In contrast, all replicashad at least one tunneling in RPM. Most of them had more than two tunnelings. Therefore,RPM samples the temperature space efficiently than REM.It was reported in experiments that C-peptide has a helix structure with salt bridges(SBs) between Glu2 and Arg10 at a low temperature [33, 34]. We obtained helix structureswhich had such SBs in our RPMD and REMD simulations, as in these reports. The lowestpotential-energy conformation among these helix structures for each simulation is presentedin Fig. 9. Here, we employed the DSSP (define secondary structure of proteins) criteria [48]for hydrogen bonds between the side-chains of Glu2 and Arg10 and for secondary structuresof C-peptide. The structure in Fig. 9(a) obtained from the RPMD simulation has twohydrogen bonds between O ǫ of Glu2 and H η of Arg10 and between O ǫ of Glu2 and H ǫ ofArg10. The residues from Ala4 to Ala11 form the α -helix structure. As for the structure18rom the REMD simulation in Fig. 9(b), the two O ǫ atoms of Glu2 form the hydrogenbonds with the two H η atoms of Arg10. The α -helix structure is formed between Ala4 andLeu9. Although these structures are slightly different with each other, both RPM and REMsampled conformations near the lowest potential-energy helix structure in the other method.To see effects of the SBs on the α -helix structures, we calculated probabilities of the α -helix structures with the SBs as well as without them. These probabilities at 281 K foreach residue are shown in Fig. 10. The probabilities without the SBs in RPM agree wellwith those in REM. The probabilities for residues 4 to 7 in both methods are high regardlessof the existence of the SBs. This is because their amino-acid sequence is AAAK, and thissequence is known for having α -helix structures [49]. In RPM the probabilities with the SBsare especially higher than those without the SBs while both probabilities are almost thesame in REM. We discuss the origin of this difference between RPM and REM later.It is considered that the SBs between Glu2 and Arg10 stabilize the α -helix structure ofC-peptide [34]. In order to investigate the relation between the SBs and the stability ofthe α -helix structure, we calculated a free-energy landscape at 281 K in each method fromEq. (36) with the reweighting techniques. These free-energy landscapes are shown in Fig. 11.The abscissa is the dihedral-angle distance d α with respect to a reference α -helix structure.Here, a dihedral-angle distance d α is defined by d α = 1 nπ n X i =1 δ ( v i , v i ) , (41)where n is the total number of dihedral angles, v i is the dihedral angle i , and v i is thedihedral angle i of the reference conformation. The distance δ ( v i , v i ) between two dihedralangles is given by δ ( v i , v i ) = min( | v i − v i | , π − | v i − v i | ) . (42)For d α , only the backbone-dihedral angles of the residues 4-10 were employed as the elementsin Eq. (41). We set the value of v i to ( φ, ψ ) = ( − π/ , − π/ d α is closeto 0, therefore, C-peptide has a α -helix structure. The ordinate is distance between O ǫ ofGlu2 and H η of Arg10: D (E2O ǫ − R10H η ) . Here, this distance is defined to be the smallestamong the four sets of the distances between two O ǫ atoms of Glu2 and two H η atoms ofArg10. Six local-minimum free-energy states are observed in both RPM and REM as shownin Figs. 11(a) and (b). We label these local-minimum states as A to F for RPM and A ′
19o F ′ for REM. The transition states between States A and B and between States A ′ andB ′ , are labelled as G and G ′ , respectively. The α -helix structure with the SBs as in Fig. 9corresponds to a structure in State A or A ′ . This fact means that the α -helix structurewith the SBs is a stable structure. On the other hand, the α -helix structure without theSBs is not a stable structure because there is no local-minimum states for this structure.Therefore, the SBs play an important role in stabilizing the α -helix structure.The global-minimum state at 281 K in RPM is State A while that in REM is StateB ′ . To see the reason for this difference, we counted the number of visiting times in StateA for RPM and in State A ′ for REM at the lowest temperature during the simulations.The region from 0.00 to 0.17 for d α and from 1.0 ˚A to 2.2 ˚A for D (E2O ǫ − R10H η ) wasassigned to State A and State A ′ here. The number of visiting times for each method islisted in Table V. Here, we employed two criteria to count the number of visiting times. InCriterion 1, when a replica visited in State A (or A ′ ) at the lowest temperature after thereplica had sampled d α larger than 0.25 or D (E2O ǫ − R10H η ) larger than 3.2 ˚A at the lowesttemperature, it is counted as one visit. In Criterion 2, sampling D (E2O ǫ − R10H η ) largerthan 3.2 ˚A was not taken into account. Therefore, the number of visiting times increasesonly when a conformation is changed from a non-helical structure to a helix structure withthe SBs in Criterion 2. In addition to this, breaking and forming of the SBs also increasesthe number of visiting times in Criterion 1. As shown in Table V, the numbers of visitingtimes in RPM are much larger than those in REM in both criteria. In Criterion 2, especially,the number of visiting times is 1 in REM. Because of such insufficient sampling in State A ′ ,this state was underestimated in REM. This underestimation caused the small probabilitiesof the α -helix structures with the SBs in Fig. 10(b).The representative conformations for each state in Fig. 11(a) are shown in Fig. 12. Therepresentative conformation for the transition state G is also presented. From these con-formations and the free-energy landscape, we clarify a folding pathway from an extendedstructure in State F to the α -helix structure with the SBs in State A. Because the pathwaysin RPM and REM are almost the same, we will discuss that only in RPM: (State F toState C via States E and D) The extended structure of C-peptide in State F changes to aglobular structure as the side-chains of Gly2 and Arg10 get close together. Turn structuresare formed between Gly2 and Arg10 by this conformational change, as in States D and E. InState C the antiparallel β -bridge structure is occasionally created between Gly2 and Arg10.20State C to State B) The SBs are formed between Gly2 and Arg10 by coming close together.Antiparallel β -bridges between Thr3 and Lys7 or between Leu9 and His12 are occasionallyobserved in State B. (State B to State A via State G) In State G, a short α -helix or 3 -helixstructure is formed around residues 7 to 10, while maintaining the SBs. By growing thishelix structure, the longer α -helix structure with the SBs is formed between residues 4 and11 as in State A.This folding pathway can be divided into two steps. The first step is the “salt-bridgeformation step”, and the second step is the “ α -helix formation step”. The first step andsecond step correspond to the transition from State F to B and that from B to A, respectively.These steps are drawn by the white and red arrows in Figs. 11 and 12. We can also seein Fig. 11 that C-peptide rarely takes a folding pathway in which the salt-bridge is formedafter the α -helix structure formation. V. CONCLUSIONS
We proposed the replica-permutation method (RPM), in which the replicas are allowedto transit not only neighboring temperatures but also non-neighboring temperatures. Forreplica-permutation trials in this method, the Suwa-Todo algorithm was employed instead ofthe Metropolis algorithm. This is because most of the permutation trials are rejected in theMetropolis algorithm. The Suwa-Todo algorithm had been proposed originally to minimizeaverage rejection ratios for state transitions in MCMC. We applied RPM and M-RPM to theparticles with the double-well potential energy to clarify the usefulness of the Suwa-Todoalgorithm for the replica permutations. For a comparison purposes, REMD simulations werealso performed. As a result, RPM realized the most efficient sampling in the temperaturespace while replica permutations were hardly accepted in M-RPM.We also applied RPM and REM to Met-enkephalin in vacuum. RPM sampled the tem-perature space more efficiently than REM even in the biomolecular system. The five local-minimum free-energy states were obtained at 200 K in both methods by using the reweightingtechniques. In the free-energy landscape estimated from the raw histogram in REM, how-ever, two of the five local-minimum states were not observed. This is because REM did notsample these two states at the lowest temperature. On the other hand, RPM sampled allstates even at the lowest temperature. It indicates that RPM realized efficient sampling not21nly in the temperature space but also in the conformational space.Furthermore, the RPMD and REMD simulations were performed for the C-peptide analogin explicit water to see the usability of RPM for a larger biomolecular system. RPM showedhigher sampling efficiency in the temperature space, again.It is reported in experiments that C-peptide has the α -helix structure with the SBsbetween Gly2 and Arg10. We observed the α -helix structures in both simulations. Wealso showed that the SBs play an important role in stabilizing the α -helix structure. Fromthe free-energy landscape, furthermore, the folding pathway from the extended structureto the α -helix structure with the SBs was clarified. This folding pathway consists of thetwo steps. The first step is the “salt-bridge formation step”. In this step, the SBs areformed by changing its conformation from the extended structure to the globular structures.The second step is the “ α -helix formation step”. The α -helix structure is formed whilemaintaining the SBs.We thus revealed that RPM realizes more efficient sampling in the conformational space atthe low temperature than REM. Furthermore, because the transition ratios of the replicas inRPM were larger than those in REM at all temperatures for all systems, larger temperatureintervals can be taken in RPM. Therefore, the number of replicas can be reduced.Although only the results of the MD simulations of RPM were shown in this article, thismethod can be readily applied to the MC method. It is also straightforward to introducethis method to the multidimensional REM [50] (also called Hamiltonian REM [51]) andrelated methods [52]. We can enhance the sampling efficiency of these methods by replacingREM to RPM. ACKNOWLEDGMENTS
The computations were performed on the computers at the Research Center for Compu-tational Science, Okazaki Research Facilities, National Institutes of Natural Sciences. [1] Mitsutake, A.; Sugita, Y.; Okamoto, Y.
Biopolymers , , 96–123.[2] Itoh, S. G.; Okumura, H.; Okamoto, Y. Mol. Sim. , , 47–56.
3] Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E.
J. Chem.Phys. , , 1087–1092.[4] Hoover, W. G.; Ladd, A. J. C.;Moran, B. Phys. Rev. Lett. , , 1818–1820.[5] Evans, D. J. J. Chem. Phys. , , 3297-3302.[6] Nos´e, S. Mol. Phys. , , 255–268.[7] Nos´e, S. J. Chem. Phys. , , 511–519.[8] Hoover, W. G. Phys. Rev. A , , 1695–1697.[9] Okumura, H.; Itoh, S. G.; Okamoto, Y. J. Chem. Phys. , , 084103.[10] Hukushima, K.; Nemoto, K. J. Phys. Soc. Jpn. , , 1604–1608.[11] Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. , , 141–151.[12] Hansmann, U. H. E.; Okamoto, Y.; Eisenmenger, F. Chem. Phys. Lett. . , 321–330.[13] Nakajima, N.; Nakamura, H.; Kidera, A. J. Phys. Chem. B , ⁀ Phys. Lett. B , , 249–253.[15] Berg, B. A.; Neuhaus, T. Phys. Rev. Lett. , , 9–12.[16] Okumura, H.; Okamoto, Y. Chem. Phys. Lett. , , 391–396.[17] Okumura, H.; Okamoto, Y. Phys. Rev. E , , 026702.[18] Okumura, H.; Okamoto, Y. J. Phys. Soc. Jpn. , , 3304–3311.[19] Okumura, H.; Okamoto, Y. Chem. Phys. Lett. , , 248–253.[20] Okumura, H.; Okamoto, Y. J. Comput. Chem. , , 379–395.[21] Berg, B. A.; Noguchi, H.; Okamoto, Y. Phys. Rev. E , , 036126.[22] Itoh, S. G.; Okamoto, Y. Mol. Sim. , , 83–89.[23] Okumura, H. J. Chem. Phys. , , 124116.[24] Okumura, H. Phys. Chem. Chem. Phys. , , 114–126.[25] Okumura, H. Proteins , in press.[26] Suwa, H.; Todo, S.
Phys. Rev. Lett. , , 120603.[27] Mitsutake, A.; Hansmann, U.H.E.; Okamoto, Y. J. Mol. Graph. Model. , , 226–238.[28] Itoh, S. G.; Okamoto, Y. Chem. Phys. Lett. , , 308–313.[29] Itoh, S. G.; Okamoto, Y. J. Chem. Phys. , , 104103.[30] Itoh, S. G.; Okamoto, Y. Phys. Rev. E , , 026705.[31] Shoemaker, K. R.; Kim, P. S.; Brems, D. N.; Marqusee, S.; York, E. J.; Chaiken, I. M.;Stewart, J. M.; Baldwin, R. L. Proc. Natl. Acad. Sci. USA , , 2349–2353.
32] Shoemaker, K. R.; Kim, P. S.; York, E. J.; Stewart, J. M.; Baldwin, R. L.
Nature , ,563–567.[33] Osterhout Jr., J. J.; Baldwin, R. L.; York, E. J.; Stewart, J. M.; Dyson, H. J.; Wright, P. E. Biochemistry , , 7059–7064.[34] Fairman, R.; Shoemaker, K. R.; York, E. J.; Stewart, J. M.; Baldwin R. L. Biophys. Chem. , , 107–119.[35] Lim, D.; Moye-Sherman, D.; Ham, I.; Jin, S.; Scholtz, J. M.; Burgess K. Chem. Commun. , 2375–2376.[36] Ferrenberg, A. M.; Swendsen, R. H.
Phys. Rev. Lett. , , 1195–1198.[37] Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput.Chem. , , , 1011–1021.[38] Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz Jr., K. M.; Ferguson, D. M.;Spellmeyer, D. C.; Fox, T.; Caldwell, J. W,; Kollman, P. A. J. Am. Chem. Soc. , ,5179–5197.[39] Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Proteins , , 712–725.[40] van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Phys. , , 1311–1327.[41] Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem.Phys. , , 926–935.[42] Ewald, P. Ann. Phys. , , 253–287.[43] Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press;Oxford, 1987.[44] MacKerell, Jr., A. D.; Bashford, D.; Bellott, M.; Dunbrack, Jr., R. L.; Evanseck, J. D.; Field,M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera,K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher,III, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.;Wi´orkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B , , 3586–3616.[45] Quenouille, M. H. Biometrika , , 353–360.[46] Miller, R. G. Biometrika , , 1–15.[47] Berg, B. A. Markov Chain Monte Carlo Simulations and Their Statistical Analysis; WorldScientific; Singapore, 2004.
48] Kabsch, W.; Sander, C.
Biopolymers , , 2577–2637.[49] Acevedo, O. E.; Lareo, L. R. OMICS: J. Integrat. Biol. , , 391–399.[50] Sugita, Y.; Kitao, A.; Okamoto, Y. J. Chem. Phys. , , 6042–6051.[51] Fukunishi, H.; Watanabe, O.; Takada, S. J. Chem. Phys. , , 9058–9067.[52] Itoh, S. G.; Okumura, H.; Okamoto, Y. J. Chem. Phys. , , 134105.[53] Sayle, R. A.; Milner-White, E. J. Trends Biochem. Sci. , ⁀
20, 374–376. ABLE I: Example of an assignment of the labels in RPM with 3 replicas. ❳❳❳❳❳❳❳❳❳❳❳❳❳❳❳❳❳❳
Temperature Label α of X α T Replica 1 Replica 1 Replica 2 Replica 2 Replica 3 Replica 3 T Replica 2 Replica 3 Replica 1 Replica 3 Replica 1 Replica 2 T Replica 3 Replica 2 Replica 3 Replica 1 Replica 2 Replica 1TABLE II: The total numbers of tunneling times during the simulations.Double-Well ∗ Enkephalin C-peptideRPM 424 ±
64 459 ±
21 58 ± ±
39 357 ±
14 27 ± ± ∗ The values were obtained by taking an average of the 40 simulations’ resultsTABLE III: The number of visiting times in each state of Met-enkephalin in vacuum.Method A B C D ERPM 97 ± ± ± ± ± ± ± ± ± ± ABLE IV: RMSD ranges for each local-minimum free-energy state.State RMSD (˚A) RMSD (˚A)A 0.0 - 0.6 1.3 - 1.8B 1.3 - 1.7 0.4 - 1.0C 1.5 - 1.8 1.3 - 1.7D 2.1 - 2.4 1.8 - 2.1E 2.0 - 2.4 2.3 - 2.8TABLE V: The numbers of visiting times in State A for RPM and in State A ′ for REM.Method Criterion 1 Criterion 2RPM 14 ± ± ± ± (cid:20)(cid:3) (cid:90) (cid:21)(cid:3) (cid:90) (cid:22)(cid:3) (cid:90) (cid:23)(cid:3) (cid:90) (cid:20)(cid:3) (cid:90) (cid:21)(cid:3) (cid:90) (cid:22)(cid:3) (cid:90) (cid:23)(cid:3) (cid:90) (cid:20)(cid:3) (cid:90) (cid:21)(cid:3) (cid:90) (cid:22)(cid:3) (cid:90) (cid:23)(cid:3) (cid:89)(cid:11)(cid:20)(cid:3450)(cid:21)(cid:12)(cid:3)(cid:89)(cid:11)(cid:20)(cid:3450)(cid:22)(cid:12)(cid:3)(cid:89)(cid:11)(cid:21)(cid:3450)(cid:22)(cid:12)(cid:3)(cid:89)(cid:11)(cid:21)(cid:3450)(cid:23)(cid:12)(cid:3)(cid:89)(cid:11)(cid:22)(cid:3450)(cid:20)(cid:12)(cid:3)(cid:89)(cid:11)(cid:22)(cid:3450)(cid:23)(cid:12)(cid:3)(cid:89)(cid:11)(cid:23)(cid:3450)(cid:20)(cid:12)(cid:3) (cid:89)(cid:11)(cid:20)(cid:3450)(cid:21)(cid:12)(cid:3)(cid:89)(cid:11)(cid:20)(cid:3450)(cid:22)(cid:12)(cid:3)(cid:89)(cid:11)(cid:21)(cid:3450)(cid:22)(cid:12)(cid:3)(cid:89)(cid:11)(cid:21)(cid:3450)(cid:23)(cid:12)(cid:3)(cid:89)(cid:11)(cid:22)(cid:3450)(cid:20)(cid:12)(cid:3) (cid:89)(cid:11)(cid:22)(cid:3450)(cid:23)(cid:12)(cid:3)(cid:89)(cid:11)(cid:23)(cid:3450)(cid:20)(cid:12)(cid:3) (cid:90) (cid:20)(cid:3) (cid:90) (cid:21)(cid:3) (cid:90) (cid:22)(cid:3) (cid:90) (cid:23)(cid:3) (cid:89)(cid:11)(cid:20)(cid:3450)(cid:20)(cid:12) (cid:3) (cid:89)(cid:11)(cid:21)(cid:3450)(cid:20)(cid:12)(cid:3) (cid:89)(cid:11)(cid:23)(cid:3450)(cid:20)(cid:12)(cid:3) (cid:89)(cid:11)(cid:20)(cid:3450)(cid:21)(cid:12) (cid:3) (cid:89)(cid:11)(cid:20)(cid:3450)(cid:22)(cid:12) (cid:3) (cid:89)(cid:11)(cid:20)(cid:3450)(cid:23)(cid:12)(cid:3)(cid:89)(cid:11)(cid:21)(cid:3450)(cid:22)(cid:12)(cid:3)(cid:89)(cid:11)(cid:21)(cid:3450)(cid:23)(cid:12)(cid:3)(cid:89)(cid:11)(cid:22)(cid:3450)(cid:20)(cid:12)(cid:3)(cid:89)(cid:11)(cid:22)(cid:3450)(cid:23)(cid:12)(cid:3) (cid:89)(cid:11)(cid:22)(cid:3450)(cid:22)(cid:12)(cid:3) (cid:89)(cid:11)(cid:22)(cid:3450)(cid:21)(cid:12)(cid:3)(cid:89)(cid:11)(cid:23)(cid:3450)(cid:21)(cid:12)(cid:3) (cid:89)(cid:11)(cid:23)(cid:3450)(cid:22)(cid:12)(cid:3) (cid:89)(cid:11)(cid:23)(cid:3450)(cid:23)(cid:12)(cid:3) (cid:90) (cid:20)(cid:3) (cid:90) (cid:21)(cid:3) (cid:90) (cid:22)(cid:3) (cid:90) (cid:23)(cid:3) (cid:89)(cid:11)(cid:21)(cid:3450)(cid:20)(cid:12)(cid:3) (cid:89)(cid:11)(cid:23)(cid:3450)(cid:20)(cid:12)(cid:3)(cid:89)(cid:11)(cid:20)(cid:3450)(cid:21)(cid:12)(cid:3)(cid:89)(cid:11)(cid:20)(cid:3450)(cid:22)(cid:12)(cid:3)(cid:89)(cid:11)(cid:20)(cid:3450)(cid:23)(cid:12)(cid:3)(cid:89)(cid:11)(cid:21)(cid:3450)(cid:22)(cid:12)(cid:3)(cid:89)(cid:11)(cid:21)(cid:3450)(cid:23)(cid:12)(cid:3)(cid:89)(cid:11)(cid:22)(cid:3450)(cid:20)(cid:12)(cid:3) (cid:89)(cid:11)(cid:22)(cid:3450)(cid:23)(cid:12)(cid:3)(cid:89)(cid:11)(cid:22)(cid:3450)(cid:22)(cid:12)(cid:3) (cid:89)(cid:11)(cid:22)(cid:3450)(cid:21)(cid:12)(cid:3)(cid:89)(cid:11)(cid:23)(cid:3450)(cid:21)(cid:12)(cid:3)(cid:89)(cid:11)(cid:23)(cid:3450)(cid:22)(cid:12)(cid:3)(cid:89)(cid:11)(cid:23)(cid:3450)(cid:23)(cid:12)(cid:3) (cid:90) (cid:20)(cid:3) (cid:90) (cid:21)(cid:3) (cid:90) (cid:22)(cid:3) (cid:90) (cid:23)(cid:3) (cid:89)(cid:11)(cid:20)(cid:3450)(cid:20)(cid:12)(cid:3) (a)(b) FIG. 1: Schematic figures of the weight allocation of the (a) Metropolis and (b) Suwa-Todoalgorithms. Red frame blocks represent rejected flows v ( i → i ) ( i = 1 , · · · , T T T T Replica 3 Replica 1
Replica 2
Replica 4 Time T e m p er a t u re Only in RPM
FIG. 2: An example of time series of temperatures in RPM. The transitions of replicas in the redsquare frame is not realized in REM. .000.030.060.090.12−3 −2 −1 0 1 2 3 P ( q ) q (a) RPMT = 200 K −2 −1 0 1 2 3q (b) RPMT = 450 K P ( q ) q (c) REMT = 200 K −2 −1 0 1 2 3q (d) REMT = 450 K P ( q ) q (e) M−RPMT = 200 K −2 −1 0 1 2 3q (f) M−RPMT = 450 K FIG. 3: The probability distributions P ( q ) of the coordinates q at (a), (c), (e) T = 200 K and(b), (d), (f) T = 450 K. These results were obtained from the (a), (b) RPMD simulations, (c),(d) REMD simulations, and (e), (f) M-RPMD simulations. The solid lines are the probabilitydistributions in Eq. (38). T r an s i t i on R a t i o Temperature Index (a) Double−Well
RPMREMM−RPM 0.00.20.40.60.8 1 2 4 6 8 10 12 T r an s i t i on R a t i o Temperature Index (b) Met−enkephalin
RPMREM 0.00.20.40.6 1 4 8 12 16 20 24 T r an s i t i on R a t i o Temperature Index (c) C−peptide
RPMREM
FIG. 4: Transition ratios of the replicas for (a) particles in the double-well potential, (b) Met-enkephalin in vacuum, and (c) C-peptide in explicit water. Temperatures are represented as thetemperature indices. The smallest and the highest indices correspond the lowest and the highesttemperatures, respectively. ∆ q ( t ) ( − ) t (ps) (a) T = 200 K M−RPMREMRPM ∆ q ( t ) ( − ) t (ps) (b) T = 450 K M−RPMREMRPM
FIG. 5: Average deviation ∆ q ( t ) in Eq. (39) at (a) T = 200 K and (b) T = 450 K. The solid line,the dashed line, and the dotted line show ∆ q ( t ) obtained from the RPMD simulations, the REMDsimulations, and the M-RPMD simulations, respectively. (a) (b) NC NC
FIG. 6: The (a) global-minimum and (b) local-minimum structures of Met-enkephalin in vacuumfor the CHARMM22 force field. MSD ( Å ) R M S D ( Å ) (a)Reweighted
0 1 2 3 0 1 2 3 A BC DE 0 2 4 6 8
RMSD ( Å ) R M S D ( Å ) (b) Reweighted
0 1 2 3 0 1 2 3 A BC DE
RMSD ( Å ) R M S D ( Å ) (c) Raw
0 1 2 3 0 1 2 3 A BC DE 0 2 4 6 8
RMSD ( Å ) R M S D ( Å ) (d) Raw
0 1 2 3 0 1 2 3 A BC DEN
A B CD E
NC NC NC NC C
FIG. 7: Free-energy landscapes at T = 200 K obtained by the reweighting techniques from the (a)RPMD and (b) REMD simulations and those calculated from the raw histograms obtained by the(c) RPMD and (d) REMD simulations. The abscissa is the RMSD with respect to the structurein Fig 6(a). The ordinate is the RMSD with respect to the structure in Fig 6(b). The unit of thefree-energy landscape is kcal/mol. The labels A to E show the global-minimum and local-minimumfree-energy states. The representative conformations at these states are also presented. The dottedlines denote hydrogen bonds. The figures were drawn by RasMol [53]. T e m pe r a t u r e I nde x Time (ns)
RPM
Replica 1 T e m pe r a t u r e I nde x Time (ns)
REM
Replica 1 T e m pe r a t u r e I nde x Time (ns)Replica 9 T e m pe r a t u r e I nde x Time (ns)Replica 9 T e m pe r a t u r e I nde x Time (ns)Replica 17 T e m pe r a t u r e I nde x Time (ns)Replica 17 FIG. 8: Time series of the temperatures of Replica 1, Replica 9, and Replica 17. The left-hand figures and the right-hand figures are obtained from the RPMD simulation and the REMDsimulation, respectively. The temperatures are represented as the temperature indices. Index 1and 24 correspond the lowest and highest temperatures, respectively. The production runs startedfrom the green dashed lines. Red circles indicate the steps at which the replicas reached the highest(lowest) temperature after they had visited the lowest (highest) temperature. The red numberspresent the tunneling times of the replicas. a) (b) NC NC
FIG. 9: The lowest potential-energy conformations of C-peptide, which had α -helix structureswith SBs between Glu2 and Arg10, obtained from the (a) RPMD and (b) REMD simulations. Thedotted lines denote the hydrogen bonds between the side-chains. The figures were created withRasMol [53]. P r obab ili t y Residue
RPM(a)
SBno SB 0.00.10.20.3 1 5 9 13 P r obab ili t y Residue
REM(b)
SBno SB
FIG. 10: Probabilities of α -helix structures at 281 K for each residue in the (a) RPMD and (b)REMD simulations. Solid line and dashed line show the results with and without SBs betweenGlu2 and Arg10, respectively. a) RPM d α D ( E O ε − R H η ) ( Å ) A BCDE FG (b) REM d α D ( E O ε − R H η ) ( Å ) A’ B’C’D’E’F’G’
FIG. 11: Free-energy landscapes at 281 K obtained from the (a) RPMD and (b) REMD simula-tions. The abscissa is the dihedral-angle distance with respect to the reference α -helix structure.The ordinate is distance between O ǫ of Glu2 and H η of Arg10. The unit of the free-energy land-scapes is kcal/mol. The labels A to F and A ′ to F ′ show the local-minimum free-energy states.The labels G and G ′ are the transition states between States A and B and between States A ′ andB ′ , respectively. A BDE