text:
# 1 Introduction

## 1 Introduction

The vertex representing a gluon splitting into a heavy-quark pair, $`g`$ $``$$`Q\overline{Q}`$ ($`Q`$ = $`b`$ or $`c`$), is a fundamental elementary component of Quantum Chromodynamics (QCD), but the contribution of this vertex to physical processes is poorly known, both theoretically and experimentally. In high-energy e+e-annihilation the leading-order process containing this vertex is e+e-$``$$`q\overline{q}g`$ $``$$`q\overline{q}`$ $`Q\overline{Q}`$ ($`q`$ = $`u,d,s,c`$ or $`b`$). Information on $`g`$ $``$$`Q\overline{Q}`$ can hence be obtained by studying e+e-$``$hadrons events comprising four quark  jets, with two of the jets identified as $`Q`$ or $`\overline{Q}`$. Background events of the kind e+e-$``$$`Q\overline{Q}`$ $`g`$ $``$$`Q\overline{Q}`$ $`q\overline{q}`$ are in principle indistinguishable final states, although their kinematics are typically quite different: the $`Q`$ and $`\overline{Q}`$ jets in this process tend to be back-to-back, and have energy comparable with the beam energy. By contrast, the $`Q`$ and $`\overline{Q}`$ jets from $`g`$ $``$$`Q\overline{Q}`$ tend to be collinear, and have low energy. In order to enhance the contribution from the latter process it is necessary to identify four-jet events and to tag typically the two jets closest in angle, and/or of lowest energy, as $`Q`$ or $`\overline{Q}`$.

We define the rate $`g_{Q\overline{Q}}`$ as the fraction of e+e-$``$hadrons events in which a gluon splits into $`Q\overline{Q}`$ , e+e-$``$$`q\overline{q}g`$ $``$$`q\overline{q}`$ $`Q\overline{Q}`$ . Since the quark mass provides a natural cutoff $`g_{Q\overline{Q}}`$ is an infrared finite quantity, which can be computed in the framework of perturbative QCD. At the $`Z^0`$ resonance energy the production of secondary $`c\overline{c}`$ or $`b\overline{b}`$ via gluon splitting is strongly suppressed by the required large gluon virtual mass. From a leading-order + next-to-leading-logarithm approximation calculation one expects  $`g_{c\overline{c}}`$ to be at the per cent level, and $`g_{b\overline{b}}`$ to be at the 0.1 per cent level. The calculations, however, depend on $`\alpha _s`$and on the quark mass, which results in substantial theoretical uncertainties.

Here we consider the process $`g`$ $``$$`b\overline{b}`$ . The measurement of $`g_{b\overline{b}}`$ is difficult experimentally since the rate is intrinsically low and the backgrounds from $`Z^0`$ $``$$`b\overline{b}`$ events are two orders of magnitude larger. In addition, the $`B`$ hadrons from $`g`$ $``$$`b\overline{b}`$ have relatively low energy and short flight distance and are difficult to identify using standard tagging techniques. So far, measurements of $`g_{b\overline{b}}`$ have been reported by ALEPH, DELPHI and OPAL ; the most precise measurement, from ALEPH, is $`g_{b\overline{b}}`$ = $`\left(2.77\pm 0.71\right)\times 10^3`$.

Such limited knowledge of $`g_{b\overline{b}}`$ results in the main source of uncertainty in the measurement of the partial decay width $`R_b=\mathrm{\Gamma }\left(Z^0b\overline{b}\right)/\mathrm{\Gamma }\left(Z^0q\overline{q}\right)`$ , which is potentially sensitive, via loop effects, to new physics processes that couple to the $`b`$-quark. Hence, more precise measurements of $`g_{b\overline{b}}`$ would help to improve the precision of tests of the electroweak theory in the heavy-quark sector. In addition, knowledge of the $`g`$ $``$$`b\overline{b}`$ process is vital for measurements at hadron colliders. For example, about $`50`$% of the $`B`$ hadrons produced in QCD processes at the Tevatron are due to $`g`$ $``$$`b\overline{b}`$ , and a larger fraction is expected to contribute at the LHC. These events form a large background to possible rare new processes involving decays to heavy quarks, such as $`H^0`$ $``$$`b\overline{b}`$ , and improved understanding of $`g`$ $``$$`b\overline{b}`$ will help to constrain background heavy-flavor production at hadron colliders.

We present a measurement of $`g_{b\overline{b}}`$ based on the sample of roughly 400,000 hadronic $`Z^0`$ decays produced in e+e-annihilations at the Stanford Linear Collider (SLC) between 1996 and 1998 and collected in the SLC Large Detector (SLD). In this period, $`Z^0`$ decays were collected with an upgraded vertex detector with wide acceptance and excellent impact parameter resolution, thus improving considerably our tagging capability for the low-energy $`B`$ hadrons characteristic of the $`g`$ $``$$`b\overline{b}`$ process.

## 2 The SLD

A description of the SLD is given elsewhere . Only the details most relevant to this analysis are mentioned here. The trigger and selection criteria for $`Z^0`$ $``$hadrons events are described elsewhere . This analysis used charged tracks measured in the Central Drift Chamber (CDC)  and in the upgraded CCD Vertex Detector (VXD) , with a momentum resolution of $`\sigma _p_{}/p_{}`$ = $`0.010.0026p_{}`$, where $`p_{}`$ is the track transverse momentum with respect to the beamline, in GeV/$`c`$. For high-momentum tracks the measured impact-parameter resolution approaches $`7.7\mu `$m ($`9.6\mu `$m) in the plane transverse to (containing) the beamline, while multiple scattering contributions are $`29\mu \mathrm{m}/\left(p\mathrm{sin}^{3/2}\theta \right)`$ in both projections. In e+e-$``$hadrons events the centroid of the SLC interaction point (IP) was reconstructed with a precision of approximately 5$`\mu `$m (10$`\mu `$m).

Only well-reconstructed tracks  were used for $`B`$-hadron tagging, and tracks from identified $`\gamma `$ conversions and $`K^0`$ or $`\mathrm{\Lambda }^0`$ decays were removed from consideration. Each track was required to have: a polar angle satisfying $`\left|\mathrm{cos}\theta \right|`$ $`<`$ 0.87, a transverse impact parameter $`<`$ 0.30cm with an error $`<`$ 250$`\mu `$m, impact parameters, measured in the CDC only, of $`<`$ 1.0cm (transverse) and $`<`$ 1.5cm (plane containing the beamline), at least 23 hits in the CDC with the first hit $`<`$ 50.0cm from the IP, a $`\chi ^2/d.o.f.`$ $`<`$ 8.0 for the CDC-only and CDC+VXD fits, and $`p_{}`$ $`>`$ 0.25GeV/c.

For the purpose of estimating the efficiency and purity of the $`g`$ $``$$`b\overline{b}`$ selection procedure, we made use of a detailed Monte-Carlo simulation of the detector. The JETSET 7.4  event generator was used, with parameter values tuned to hadronic $`e^+e^{}`$ annihilation data , combined with a simulation of $`B`$ hadron decays tuned to $`\mathrm{{\rm Y}}\left(4S\right)`$ data  and a simulation of the SLD based on GEANT 3.21 . Inclusive distributions of single-particle and event-topology observables in hadronic events were found to be well described by the simulations . Uncertainties in the simulation were taken into account in the systematic errors (Section 6).

Monte-Carlo events were reweighted to take into account the current measurements of gluon splitting into heavy-quark pairs . JETSET with the SLD parameters predicts $`g_{b\overline{b}}=0.14\%`$ and $`g_{c\overline{c}}=1.36\%`$. We reweighted these events in the simulated sample to obtain $`g_{b\overline{b}}=0.247\%`$ and $`g_{c\overline{c}}=3.07\%`$ . Samples of about 1900k Monte-Carlo $`Z^0`$ $``$$`q\overline{q}`$ events, 1900k $`Z^0`$ $``$$`b\overline{b}`$ events, 1090k $`Z^0`$ $``$$`c\overline{c}`$ events and 60k $`Z^0`$ $``$$`q\overline{q}g`$ , $`g`$ $``$$`b\overline{b}`$ events were used to evaluate the selection efficiencies (Section 6).

## 3 Flavor Tagging

We used topologically-reconstructed secondary vertices  for heavy-quark tagging. To reconstruct the secondary vertices, the space points where track density functions overlap were found in three dimensions. Only the vertices that are significantly displaced from the IP were considered to be possible $`B`$\- or $`D`$-hadron decay vertices. The mass of the secondary vertex was calculated using the tracks that were associated with the vertex. We corrected the reconstructed mass to account for neutral decay products and tracks missed from the vertex. By using kinematic information from the vertex flight path and the momentum sum of the tracks associated with the secondary vertex, we calculated the $`P_T`$-corrected mass, $`M_{P_T}`$, by adding a component of missing momentum to the invariant mass, as follows:

$$M_{P_T}=\sqrt{M_{vtx}^{}{}_{}{}^{2}+P_T^2}+\left|P_T\right|.$$

where $`M_{vtx}`$ is the invariant mass of the tracks associated with the reconstructed secondary vertex and $`P_T`$ is the total transverse momentum of the vertex-associated tracks with respect to the vertex axis, which we estimated independently of the track momenta by the vector along the line joining the IP to the reconstructed vertex position. In this correction, vertexing resolution as well as the IP resolution are crucial.

With these features, topological vertex finding gives excellent $`b`$-tagging efficiency and purity. In particular, the efficiency is good even at low $`B`$-hadron energies, which is especially important for detecting $`g`$ $``$$`b\overline{b}`$ . For the selected 4-jet event sample (Section 4) we used our simulation to estimate that our mean tagging efficiency for $`b`$-jets in the $`g`$ $``$$`b\overline{b}`$ process is 67%.

## 4 Analysis

Besides the signal events which contain $`g`$ $``$$`b\overline{b}`$ , hereafter called ‘B events’, background events can be divided into two categories: 1) events in which a gluon splits to a charm quark pair, called ‘C events’, and 2) events which do not contain any gluon splitting into heavy quarks at all, hereafter called ‘Q events’. In $`g`$ $``$$`b\overline{b}`$ events the two $`B`$ hadrons from the gluon splitting tend to be produced with low energy in a collinear configuration, which allows one to discriminate the signal from background. We first required that each event contain 4 jets, and that the two jets closest in angle were tagged as $`b`$ jets on the basis that each contained a secondary vertex. We then examined additional kinematic quantities and used a neural network technique to improve the signal/background ratio.

In each event jets were formed by applying the Durham jet-finding algorithm  with $`y_{cut}=0.005`$ to the set of charged tracks; this $`y_{cut}`$ value was chosen to minimize the sum of the statistical and systematic errors on $`g_{b\overline{b}}`$. Events containing four or more jets were retained. The jet-finder was re-run on the $`>4`$ jet events with successively larger $`y_{cut}`$ values until exactly four jets were reconstructed. With this definition the 4-jet rate in the data was $`\left(14.58\pm 0.07\right)\%`$, where the error is statistical only. In the Monte-Carlo simulation the 4-jet rate was $`\left(14.47\pm 0.02\pm 0.16\right)\%`$ where the first error is statistical and the second is due to the uncertainty in the simulation of heavy-quark physics (Section 6). For the B, C and Q events the 4-jet rates predicted by the simulation are about $`60\%`$, $`38\%`$ and $`14\%`$, respectively. Each jet energy was calculated using its associated-track momenta and assuming all tracks to have the charged pion mass.

In each selected event the two jets closest in angle were considered as candidates for originating from the gluon splitting process $`g`$ $``$$`b\overline{b}`$ , and the topological vertex method was applied to them. We required both jets to contain a secondary vertex with a 3D decay length greater than $`300\mu `$m. No tag was applied to the other two jets. 1514 events were selected. In each event the tagged jets were labeled ‘1’ and ‘2’, where jet 1 contained the vertex with the greater $`M_{P_T}`$ value, $`M_{P_T1}`$, and jet 2 that with the lesser $`M_{P_T}`$ value, $`M_{P_T2}`$. The other two jets in the event were labeled ‘3’ and ‘4’, where jet 3 was more energetic than jet 4. With these requirements the selection efficiency for $`g`$ $``$$`b\overline{b}`$ events was estimated to be $`16.2\%`$, with a signal/background ratio in the selected sample of approximately 1/10. $`75\%`$ of the background came from $`Z^0`$ $``$$`b\overline{b}`$ events, $`9\%`$ from $`Z^0`$ $``$$`q\overline{q}`$ ($`qb`$) events, and the remaining $`16\%`$ from $`g`$ $``$$`c\overline{c}`$ events. In order to improve the signal/background ratio we used a neural network technique. We chose the following 9 observables as inputs to the neural network; each observable was scaled to correspond to a range between 0 and 1.

1. $`M_{P_T1}`$: $`b`$ jets typically have higher values of this quantity than $`c`$ or $`uds`$ jets. The distribution of $`M_{P_T1}`$ is shown in Figure 1.
2. $`M_{P_T2}`$: This observable has similar discriminating power. The distribution of $`M_{P_T2}`$ is shown in Figure 2.
3. $`15M_{P_T1}P_{vtx1}`$: where $`P_{vtx}`$ is the vertex momentum. This observable tends to be large for $`b`$ jets since $`B`$ decay vertices typically have higher mass than those from charm decays, and vertices resulting from $`B`$ $``$$`D`$ cascade decays have a lower momentum than those from primary $`D`$ hadrons. The distribution of $`15M_{P_T1}P_{vtx1}`$ is shown in Figure 3.
4. $`15M_{P_T2}P_{vtx2}`$: This observable also has discriminating power between signal and background events. The distribution of $`15M_{P_T2}P_{vtx2}`$ is shown in Figure 4.
5. The angle $`\theta _{12}`$ between the vertex axes of jets 1 and 2. The two jets from $`g`$ $``$$`b\overline{b}`$ tend to have $`\theta _{12}`$ $`0`$. However, in some $`Z^0`$ $``$$`b\overline{b}`$ background events a single $`b`$ jet may be split into two by the jet-finder; in these cases the two reconstructed vertices tend to have $`\mathrm{cos}\theta _{12}0.98`$. The distribution of cos$`\theta _{12}`$ is shown in Figure 5.
6. The angle $`\theta _{34}`$ between the axes of jets 3 and 4. In events containing $`g`$ $``$$`b\overline{b}`$ this tends to be near $`\pi `$, while background events tend to populate the smaller-angle region. The distribution of cos$`\theta _{34}`$ is shown in Figure 6.
7. The energy sum of jets 1 and 2, $`E_1+E_2`$: The two jets arising from $`g`$ $``$$`b\overline{b}`$ tend to have lower energy than the other two jets in the event. The distribution of $`E_1+E_2`$ is shown in Figure 7.
8. The energy sum of jets 3 and 4, $`E_3+E_4`$: This tends to be larger in signal events than in background events. The distribution of $`E_3+E_4`$ is shown in Figure 8.
9. The angle $`\alpha _{1234}`$ between the plane $`\mathrm{\Pi }_{12}`$ formed by jets 1 and 2 and the plane $`\mathrm{\Pi }_{34}`$ formed by jets 3 and 4: This variable is similar to the Bengtsson-Zerwas angle , and is useful to separate $`g`$ $``$$`b\overline{b}`$ events because the radiated virtual gluon in the process $`Z^0`$ $``$$`q\overline{q}g`$ is polarized in the three-parton plane, and this is reflected in its subsequent splitting, by favoring $`g`$ $``$$`q\overline{q}`$ emission out of this plane, i.e.$`\alpha _{1234}\pi /2`$. The distribution of cos$`\alpha _{1234}`$  is shown in Figure 9.

The measured and simulated distributions agree well for these input observables. We trained the neural network using Monte-Carlo samples of about 1800k $`Z^0`$ $``$$`q\overline{q}`$ events, 1200k $`Z^0`$ $``$$`b\overline{b}`$ events, 780k $`Z^0`$ $``$$`c\overline{c}`$ events and 50k events containing $`g`$ $``$$`b\overline{b}`$ . These samples were independent of the ones used for the selection efficiency and background studies. Figure 10 shows the distribution of the neural network output variable, $`Y`$. We retained events with $`Y`$ greater than 0.7. This value was found to minimise the total error on the final $`g_{b\overline{b}}`$ result.

## 5 Result

79 events were selected in the data. The number of background events was estimated, using the Monte Carlo simulation, to be $`41.9`$, where $`35\%`$ of the background comes from $`g`$ $``$$`c\overline{c}`$ events, $`63\%`$ from $`Z^0`$ $``$$`b\overline{b}`$ events, and the remaining $`2\%`$ from $`Z^0`$ $``$$`q\overline{q}`$ $`\left(qb\right)`$ events. Table 1 shows the selection efficiencies, relative to the selected hadronic-event sample, for the B, C and Q event categories. From these efficiencies and the fraction of events selected in the data, $`f_d=\left(2.73\pm 0.31\right)\times 10^4`$, the value of $`g_{b\overline{b}}`$ was determined:

$$g_{b\overline{b}}=\frac{f_d\left(1g_{c\overline{c}}\right)ϵ_Qg_{c\overline{c}}ϵ_C}{ϵ_Bϵ_Q}.$$
(1)
We obtained

$$g_{b\overline{b}}=\left(2.44\pm 0.59\right)\times 10^3,$$
(2)
where the error is statistical only.

## 6 Systematic Errors

The efficiencies for the three event categories were evaluated using the Monte-Carlo simulation. The limitations of the simulation in estimating these efficiencies lead to an uncertainty on the result. The error due to limited Monte-Carlo statistics in the efficiency evaluation was $`\mathrm{\Delta }g_{b\overline{b}}=\pm 0.12\times 10^3`$.

A large fraction of events remaining after the selection cuts contain $`B`$ and $`D`$ hadrons. The uncertainty in the knowledge of the physical processes in the simulation of heavy-flavor production and decays constitutes a source of systematic error. All the physical simulation parameters were varied within their allowed experimental ranges. In particular, the $`B`$ and $`D`$ hadron lifetimes, their production rates, and the mean $`B`$ hadron energy were varied following the recommendations of the LEP Heavy Flavour Working Group . We also varied the assumed form of the $`B`$ energy distribution, taking the difference between the optimised Bowler and Peterson functions  to assign an error. The uncertainties, which are typically small, are summarized in Table 2.

The dominant systematic error arises via the signal tagging efficiency, $`ϵ_B`$, from the dependence on modelling the kinematics of the split gluon: its energy $`E_g`$, its mass $`m_g`$ and the decay angle, $`\theta ^{}`$, of the two $`B`$ hadrons in their center-of-mass frame relative to the gluon direction. In our default Monte Carlo simulation the kinematics of the signal events are based on the JETSET parton shower, which is in good agreement with the theoretical predictions . In order to estimate the uncertainty on this assumption, we have investigated alternative models, namely the PYTHIA 6.136, HERWIG 6.1 , and ARIADNE 4.08  models, as well as the exact leading-order matrix-element calculation GRC4F . In each case, we generated, at the parton level, events containing $`g`$ $``$$`b\overline{b}`$ . For illustration the $`E_g`$ distributions are shown in Fig. 11. The efficiency function computed with JETSET was then reweighted by the ratio of the new model to JETSET initial distributions to obtain a new estimate of the average efficiency. The resulting $`g_{b\overline{b}}`$ values are shown in Table 3. We took the central value of this ensemble, $`g_{b\overline{b}}`$ = 2.44 $`\times 10^3`$, as our central result, and assigned a systematic error of $`\mathrm{\Delta }g_{b\overline{b}}=\pm 0.24\times 10^3`$ based on the full range of values. For illustration the efficiency and its uncertainty are shown in Fig. 11.

The dependence of the efficiency on the $`b`$-quark mass was also investigated at the generator level using the GRC4F Monte Carlo program. The variation of $`ϵ_B`$ was computed for $`b`$-quark masses between $`4.7`$ and $`5.3`$ GeV/c2. The resultant uncertainty is estimated to be $`\mathrm{\Delta }g_{b\overline{b}}=\pm 0.05\times 10^3`$. The measured uncertainty in the production fraction of $`g`$ $``$$`c\overline{c}`$ background events, $`\mathrm{\Delta }g_{c\overline{c}}=\pm 0.40\%`$, gives an error $`\mathrm{\Delta }g_{b\overline{b}}=\pm 0.12\times 10^3`$.

In the Monte-Carlo simulation charged tracks used in the topological vertex tag were rejected to reproduce better the distributions in the data. Uncertainties in the efficiencies due to this rejection were assessed by evaluating the Monte-Carlo efficiencies without the rejection algorithm. The difference in the $`g_{b\overline{b}}`$ result was taken as a symmetric systematic error, $`\mathrm{\Delta }g_{b\overline{b}}=\pm 0.04\times 10^3`$.

As cross checks we varied independently the value of $`y_{cut}`$ used in the 4-jet event selection and the value of the neural network output, $`Y`$, used to select the final event sample. In each case we found results consistent with those determined using the optimised value. Table 2 summarizes the different sources of systematic error on $`g_{b\overline{b}}`$. The total systematic error was estimated to be the sum in quadrature, $`0.34\times 10^3`$.

## 7 Summary

Using the excellent flavor-tagging capabilities of the SLD tracking system, and a new technique incorporating a multivariate neural network analysis, we have measured the probability for gluon splitting into $`b\overline{b}`$ in hadronic $`Z^0`$ decays. The result is

$$g_{b\overline{b}}=(2.44\pm 0.59(\mathrm{stat}.)\pm 0.34(\mathrm{syst}.))\times 10^3$$

where the first error is statistical and the second is the sum in quadrature of systematic effects. This represents the most precise determination of $`g_{b\overline{b}}`$ to date. Our result is consistent with previous measurements . It is also consistent with the theoretical expectation $`g_{b\overline{b}}`$ = $`1.75\pm 0.40\times 10^3`$ , where the central value corresponds to $`\alpha _s\left(M_Z^2\right)`$= 0.118 and a $`b`$-quark mass of 5.0 GeV/$`c^2`$. Finally, we found that the predictions of the models PYTHIA6.136 ($`g_{b\overline{b}}`$ = $`1.5\times 10^3`$) and HERWIG6.1 ($`g_{b\overline{b}}`$ = $`2.5\times 10^3`$), with their default parameter settings, are consistent with our measurement. The prediction of ARIADNE4.08 ($`g_{b\overline{b}}`$ = $`1.3\times 10^3`$) is slightly below our measurement.

## Acknowledgments

We thank Mike Seymour for helpful conversations. We thank the personnel of the SLAC accelerator department and the technical staffs of our collaborating institutions for their outstanding efforts on our behalf. This work was supported by the U.S. Department of Energy and National Science Foundation, the UK Particle Physics and Astronomy Research Council, the Istituto Nazionale di Fisica Nucleare of Italy, the Japan-US Cooperative Research Project on High Energy Physics, and the Korea Research Foundation.

∗∗ List of Authors

Koya Abe,(24) Kenji Abe,(15) T. Abe,(21) I. Adam,(21) H. Akimoto,(21) D. Aston,(21) K.G. Baird,(11) C. Baltay,(30) H.R. Band,(29) T.L. Barklow,(21) J.M. Bauer,(12) G. Bellodi,(17) R. Berger,(21) G. Blaylock,(11) J.R. Bogart,(21) G.R. Bower,(21) J.E. Brau,(16) M. Breidenbach,(21) W.M. Bugg,(23) D. Burke,(21) T.H. Burnett,(28) P.N. Burrows,(17) A. Calcaterra,(8) R. Cassell,(21) A. Chou,(21) H.O. Cohn,(23) J.A. Coller,(4) M.R. Convery,(21) V. Cook,(28) R.F. Cowan,(13) G. Crawford,(21) C.J.S. Damerell,(19) M. Daoudi,(21) N. de Groot,(2) R. de Sangro,(8) D.N. Dong,(13) M. Doser,(21) R. Dubois,(21) I. Erofeeva,(14) V. Eschenburg,(12) E. Etzion,(29) S. Fahey,(5) D. Falciai,(8) J.P. Fernandez,(26) K. Flood,(11) R. Frey,(16) E.L. Hart,(23) K. Hasuko,(24) S.S. Hertzbach,(11) M.E. Huffer,(21) X. Huynh,(21) M. Iwasaki,(16) D.J. Jackson,(19) P. Jacques,(20) J.A. Jaros,(21) Z.Y. Jiang,(21) A.S. Johnson,(21) J.R. Johnson,(29) R. Kajikawa,(15) M. Kalelkar,(20) H.J. Kang,(20) R.R. Kofler,(11) R.S. Kroeger,(12) M. Langston,(16) D.W.G. Leith,(21) V. Lia,(13) C. Lin,(11) G. Mancinelli,(20) S. Manly,(30) G. Mantovani,(18) T.W. Markiewicz,(21) T. Maruyama,(21) A.K. McKemey,(3) R. Messner,(21) K.C. Moffeit,(21) T.B. Moore,(30) M. Morii,(21) D. Muller,(21) V. Murzin,(14) S. Narita,(24) U. Nauenberg,(5) H. Neal,(30) G. Nesom,(17) N. Oishi,(15) D. Onoprienko,(23) L.S. Osborne,(13) R.S. Panvini,(27) C.H. Park,(22) I. Peruzzi,(8) M. Piccolo,(8) L. Piemontese,(7) R.J. Plano,(20) R. Prepost,(29) C.Y. Prescott,(21) B.N. Ratcliff,(21) J. Reidy,(12) P.L. Reinertsen,(26) L.S. Rochester,(21) P.C. Rowson,(21) J.J. Russell,(21) O.H. Saxton,(21) T. Schalk,(26) B.A. Schumm,(26) J. Schwiening,(21) V.V. Serbo,(21) G. Shapiro,(10) N.B. Sinev,(16) J.A. Snyder,(30) H. Staengle,(6) A. Stahl,(21) P. Stamer,(20) H. Steiner,(10) D. Su,(21) F. Suekane,(24) A. Sugiyama,(15) S. Suzuki,(15) M. Swartz,(9) F.E. Taylor,(13) J. Thom,(21) E. Torrence,(13) T. Usher,(21) J. Va’vra,(21) R. Verdier,(13) D.L. Wagner,(5) A.P. Waite,(21) S. Walston,(16) A.W. Weidemann,(23) E.R. Weiss,(28) J.S. Whitaker,(4) S.H. Williams,(21) S. Willocq,(11) R.J. Wilson,(6) W.J. Wisniewski,(21) J.L. Wittlin,(11) M. Woods,(21) T.R. Wright,(29) R.K. Yamamoto,(13) J. Yashima,(24) S.J. Yellin,(25) C.C. Young,(21) H. Yuta.(1)

(1)Aomori University, Aomori, 030 Japan, (2)University of Bristol, Bristol, United Kingdom, (3)Brunel University, Uxbridge, Middlesex, UB8 3PH United Kingdom, (4)Boston University, Boston, Massachusetts 02215, (5)University of Colorado, Boulder, Colorado 80309, (6)Colorado State University, Ft. Collins, Colorado 80523, (7)INFN Sezione di Ferrara and Universita di Ferrara, I-44100 Ferrara, Italy, (8)INFN Laboratori Nazionali di Frascati, I-00044 Frascati, Italy, (9)Johns Hopkins University, Baltimore, Maryland 21218-2686, (10)Lawrence Berkeley Laboratory, University of California, Berkeley, California 94720, (11)University of Massachusetts, Amherst, Massachusetts 01003, (12)University of Mississippi, University, Mississippi 38677, (13)Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, (14)Institute of Nuclear Physics, Moscow State University, 119899 Moscow, Russia, (15)Nagoya University, Chikusa-ku, Nagoya, 464 Japan, (16)University of Oregon, Eugene, Oregon 97403, (17)Oxford University, Oxford, OX1 3RH, United Kingdom, (18)INFN Sezione di Perugia and Universita di Perugia, I-06100 Perugia, Italy, (19)Rutherford Appleton Laboratory, Chilton, Didcot, Oxon OX11 0QX United Kingdom, (20)Rutgers University, Piscataway, New Jersey 08855, (21)Stanford Linear Accelerator Center, Stanford University, Stanford, California 94309, (22)Soongsil University, Seoul, Korea 156-743, (23)University of Tennessee, Knoxville, Tennessee 37996, (24)Tohoku University, Sendai, 980 Japan, (25)University of California at Santa Barbara, Santa Barbara, California 93106, (26)University of California at Santa Cruz, Santa Cruz, California 95064, (27)Vanderbilt University, Nashville,Tennessee 37235, (28)University of Washington, Seattle, Washington 98105, (29)University of Wisconsin, Madison,Wisconsin 53706, (30)Yale University, New Haven, Connecticut 06511.


format: text

text:
# 1 Introduction

## 1 Introduction

Over the past three decades, a number of powerful simulation algorithms have been introduced to the protein folding problem (for reviews see, e.g., Refs. ). For many years, the emphasis has been placed on how to find the global-minimum-energy conformation among a huge number of local-minimum states. For complete understanding of protein folding mechanism, however, the global knowledge of the configurational space is required, including the intermediate and denatured states of proteins. For this purpose, free-energy calculations are essential.

We have been advocating the uses of generalized-ensemble algorithms as the methods that meet the above requirements (for reviews see, e.g., Refs. ). In this method each state is weighted by a non-Boltzmann probability weight factor so that a random walk in potential energy space may be realized. The random walk allows the simulation to escape from any energy barrier and to sample much wider configurational space than by conventional methods. Monitoring the energy in a single simulation run, one can obtain not only the global-minimum-energy state but also canonical ensemble averages as functions of temperature by the single-histogram  and/or multiple-histogram  reweighting techniques (an extension of the multiple-histogram method is also referred to as weighted histogram analysis method (WHAM) ).

Three of the most well-known generalized-ensemble methods are perhaps multicanonical algorithm (MUCA) , simulated tempering (ST) , and replica-exchange method (REM) . (MUCA is also referred to as entropic sampling  and adaptive umbrella sampling . ST is also referred to as the method of expanded ensemble . REM is also referred to as parallel tempering . Details of literature about REM and related algorithms can be found in a recent review .) Since MUCA was first introduced to protein folding problem , various generalized-ensemble algorithms have been used in many applications in protein and related systems (see Ref.  and references therein). In particular, free-energy calculations in protein folding by generalized-ensemble algorithms were explored in Refs. .

REM has been drawing much attention recently because the probability weight factors are essentially known a priori, whereas they are not in most of other generalized-ensemble algorithms (and have to be determined by a tedius procedure). In REM a number of non-interacting copies (or replicas) of the original system at different temperatures are simulated independently and simultaneously by the conventional Monte Carlo (MC) or molecular dynamics (MD) method. Every few steps, pairs of replicas are exchanged with a specified transition probability.

We have worked out the details for the replica-exchange molecular dynamics algorithm  (see also Ref. ). We have also developed a multidimensional replica-exchange method (MREM)  (see also Refs. ). In MREM we showed that REM is not limited to tempering (or temperature exchange) and that we can also exchange parameters in the potential energy. Namely, pairs of replicas with different temperatures and/or different parameters of the potential energy are exchanged during the simulation. Important applications of MREM are free-energy calculations.

The umbrella sampling method  and free energy perturbation method, which is a special case of umbrella sampling, have been widely used to calculate the free energies in chemical processes  \- . Although the effectiveness of the umbrella sampling method is well known, its successful implementation requires a careful fine tuning. Various generalizations of the umbrella sampling method have thus been introduced to sample the potential energy surface more effectively. The $`\lambda `$-$`dynamics`$  \-  is such an example, where the coupling parameter $`\lambda `$ is treated as a dynamical variable. Another example is the multicanonical WHAM , which combines the umbrella sampling with multicanonical algorithm. We have developed yet another generalization of the umbrella sampling method (we refer to this method as replica-exchange umbrella sampling (REUS)), which is based on the multidimensional extension of the replica-exchange method .

REM is very effective and has already been used in many applications in protein systems (see Ref.  and references therein). However, REM also has a computational difficulty: As the number of degrees of freedom of the system increases, the required number of replicas also greatly increases, whereas only a single replica is simulated in MUCA or ST. This demands a lot of computer power for complex systems. Our solution to this problem is: Use REM for the weight factor determinations of MUCA or ST, which is much simpler than previous iterative methods of weight determinations, and then perform a long MUCA or ST production run. The first example is the replica-exchange multicanonical algorithm (REMUCA) . In REMUCA, a short replica-exchange simulation is performed, and the multicanonical weight factor is determined by WHAM . Another example of such a combination is the replica-exchange simulated tempering (REST) . In REST, a short replica-exchange simulation is performed, and the simulated tempering weight factor is determined by WHAM . We have introduced a further extension of REMUCA, which we refer to as multicanonical replica-exchange method (MUCAREM) . In MUCAREM, the multicanonical weight factor is first determined as in REMUCA, and then a replica-exchange multicanonical production simulation is performed with a small number of replicas.

In this article, we first describe the multidimensional replica-exchange method, a particular realization of which is the replica-exchange umbrella sampling . We then present the replica-exchange multicanonical algorithm . The effectiveness of these methods is tested with short peptide systems.

## 2 Methods

### 2.1 Multidimensional Replica-Exchange Method

We first review the original replica-exchange method (REM)  (see Ref.  for details of the molecular dynamics version).

We consider a system of $`N`$ atoms with their coordinate vectors and momentum vectors denoted by $`q\{𝒒_1,\mathrm{},𝒒_N\}`$ and $`p\{𝒑_1,\mathrm{},𝒑_N\}`$, respectively. The Hamiltonian $`H(q,p)`$ of the system is the sum of the kinetic energy $`K(p)`$ and the potential energy $`E(q)`$:

$$H(q,p)=K(p)+E(q).$$
(1)
In the canonical ensemble at temperature $`T`$ each state $`x(q,p)`$ with the Hamiltonian $`H(q,p)`$ is weighted by the Boltzmann factor:

$$W_T(x)=e^{\beta H(q,p)},$$
(2)
where the inverse temperature $`\beta `$ is defined by $`\beta =1/k_BT`$ ($`k_B`$ is the Boltzmann constant).

The generalized ensemble for REM consists of $`M`$ non-interacting copies (or, replicas) of the original system in the canonical ensemble at $`M`$ different temperatures $`T_m`$ ($`m=1,\mathrm{},M`$). We arrange the replicas so that there is always exactly one replica at each temperature. Then there is a one-to-one correspondence between replicas and temperatures; the label $`i`$ ($`i=1,\mathrm{},M`$) for replicas is a permutation of the label $`m`$ ($`m=1,\mathrm{},M`$) for temperatures, and vice versa:

$$\{\begin{array}{cc}\hfill i& =i(m)f(m),\hfill \\ \hfill m& =m(i)f^1(i),\hfill \end{array}$$
(3)
where $`f(m)`$ is a permutation function of $`m`$ and $`f^1(i)`$ is its inverse.

Let $`X=\{x_1^{[i(1)]},\mathrm{},x_M^{[i(M)]}\}=\{x_{m(1)}^{[1]},\mathrm{},x_{m(M)}^{[M]}\}`$ stand for a “state” in this generalized ensemble. Here, the superscript and the subscript in $`x_m^{[i]}`$ label the replica and the temperature, respectively. 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_m^{[i]}(q^{[i]},p^{[i]})_m.$$
(4)
Because the replicas are non-interacting, the weight factor for the state $`X`$ is given by the product of Boltzmann factors for each replica (or at each temperature):

$$\begin{array}{cc}\hfill W_{\mathrm{REM}}(X)& =\mathrm{exp}\left\{\underset{i=1}{\overset{M}{}}\beta _{m(i)}H(q^{[i]},p^{[i]})\right\},\hfill \\ & =\mathrm{exp}\left\{\underset{m=1}{\overset{M}{}}\beta _mH(q^{[i(m)]},p^{[i(m)]})\right\},\hfill \end{array}$$
(5)
where $`i(m)`$ and $`m(i)`$ are the permutation functions in Eq. (3).

We now consider exchanging a pair of replicas. Suppose we exchange replicas $`i`$ and $`j`$ which are at temperatures $`T_m`$ and $`T_n`$, respectively:

$$X=\{\mathrm{},x_m^{[i]},\mathrm{},x_n^{[j]},\mathrm{}\}X^{}=\{\mathrm{},x_m^{[j]},\mathrm{},x_n^{[i]},\mathrm{}\}.$$
(6)
The exchange of replicas can be written in more detail as

$$\{\begin{array}{cc}\hfill x_m^{[i]}(q^{[i]},p^{[i]})_m& x_m^{[j]}(q^{[j]},p^{[j]})_m,\hfill \\ \hfill x_n^{[j]}(q^{[j]},p^{[j]})_n& x_n^{[i]}(q^{[i]},p^{[i]})_n,\hfill \end{array}$$
(7)
where the momenta are uniformly rescaled according to

$$\{\begin{array}{cc}\hfill p^{[i]}& \sqrt{\frac{T_n}{T_m}}p^{[i]},\hfill \\ \hfill p^{[j]}& \sqrt{\frac{T_m}{T_n}}p^{[j]}.\hfill \end{array}$$
(8)
In order for this exchange process to converge towards the equilibrium distribution based on Eq. (5), it is sufficient to impose the detailed balance condition on the transition probability $`w(XX^{})`$:

$$W_{\mathrm{REM}}(X)w(XX^{})=W_{\mathrm{REM}}(X^{})w(X^{}X).$$
(9)
From Eqs. (1), (5), (8), and (9), we have

$$\frac{w(XX^{})}{w(X^{}X)}=\mathrm{exp}\left(\mathrm{\Delta }\right),$$
(10)
where

$`\mathrm{\Delta }`$ $`=`$ $`\beta _m\left(E\left(q^{[j]}\right)E\left(q^{[i]}\right)\right)\beta _n\left(E\left(q^{[j]}\right)E\left(q^{[i]}\right)\right),`$ (11)
$`=`$ $`\left(\beta _m\beta _n\right)\left(E\left(q^{[j]}\right)E\left(q^{[i]}\right)\right).`$ (12)
This can be satisfied, for instance, by the usual Metropolis criterion:

$$w(XX^{})w\left(x_m^{[i]}|x_n^{[j]}\right)=\{\begin{array}{cc}1,\hfill & \mathrm{for}\mathrm{\Delta }0,\hfill \\ \mathrm{exp}\left(\mathrm{\Delta }\right),\hfill & \mathrm{for}\mathrm{\Delta }>0.\hfill \end{array}$$
(13)
Note that because of the velocity rescaling of Eq. (8) the kinetic energy terms are cancelled out in Eqs. (11) (and (12)) and that the same criterion, Eqs. (12) and (13), which was originally derived for Monte Carlo algorithm  is recovered .

A simulation of the replica-exchange method (REM)  is then realized by alternately performing the following two steps:

1. Each replica in canonical ensemble of the fixed temperature is simulated $`simultaneously`$ and $`independently`$ for a certain MC or MD steps.
2. A pair of replicas, say $`x_m^{[i]}`$ and $`x_n^{[j]}`$, are exchanged with the probability $`w\left(x_m^{[i]}|x_n^{[j]}\right)`$ in Eq. (13).

In the present work, we employ molecular dynamics algorithm for Step 1. Note that in Step 2 we exchange only pairs of replicas corresponding to neighboring temperatures, because the acceptance ratio of the exchange decreases exponentially with the difference of the two $`\beta `$’s (see Eqs. (12) and (13)). Note also that whenever a replica exchange is accepted in Step 2, the permutation functions in Eq. (3) are updated.

The method is particularly suitable for parallel computers. Because one can minimize the amount of information exchanged among nodes, it is best to assign each replica to each node (exchanging pairs of temperature values among nodes is much faster than exchanging coordinates and momenta). This means that we keep track of the permutation function $`m(i;t)=f^1(i;t)`$ in Eq. (3) as a function of MD step $`t`$ throughout the simulation.

The major advantage of REM over other generalized-ensemble methods such as multicanonical algorithm  and simulated tempering  lies in the fact that the weight factor is a priori known (see Eq. (5)), whereas in the latter algorithms the determination of the weight factors can be very tedius and time-consuming. A random walk in “temperature space” is realized for each replica, which in turn induces a random walk in potential energy space. This alleviates the problem of getting trapped in states of energy local minima. In REM, however, the number of required replicas increases as the system size $`N`$ increases (according to $`\sqrt{N}`$) . This demands a lot of computer power for complex systems.

We now present our multidimensional extension of REM, which we refer to as multidimensional replica-exchange method (MREM). The crucial observation that led to the new algorithm is: As long as we have $`M`$ non-interacting replicas of the original system, the Hamiltonian $`H(q,p)`$ of the system does not have to be identical among the replicas and it can depend on a parameter with different parameter values for different replicas. Namely, we can write the Hamiltonian for the $`i`$-th replica at temperature $`T_m`$ as

$$H_m(q^{[i]},p^{[i]})=K(p^{[i]})+E_{\lambda _m}(q^{[i]}),$$
(14)
where the potential energy $`E_{\lambda _m}`$ depends on a parameter $`\lambda _m`$ and can be written as

$$E_{\lambda _m}(q^{[i]})=E_0(q^{[i]})+\lambda _mV(q^{[i]}).$$
(15)
This expression for the potential energy is often used in simulations. For instance, in umbrella sampling , $`E_0(q)`$ and $`V(q)`$ can be respectively taken as the original potential energy and the “biasing” potential energy with the coupling parameter $`\lambda _m`$. In simulations of spin systems, on the other hand, $`E_0(q)`$ and $`V(q)`$ (here, $`q`$ stands for spins) can be respectively considered as the zero-field term and the magnetization term coupled with the external field $`\lambda _m`$.

While replica $`i`$ and temperature $`T_m`$ are in one-to-one correspondence in the original REM, replica $`i`$ and “parameter set” $`\mathrm{\Lambda }_m(T_m,\lambda _m)`$ are in one-to-one correspondence in the new algorithm. Hence, the present algorithm can be considered as a multidimensional extension of the original replica-exchange method where the “parameter space” is one-dimensional (i.e., $`\mathrm{\Lambda }_m=T_m`$). Because the replicas are non-interacting, the weight factor for the state $`X`$ in this new generalized ensemble is again given by the product of Boltzmann factors for each replica (see Eq. (5)):

$$\begin{array}{cc}\hfill W_{\mathrm{MREM}}(X)& =\mathrm{exp}\left\{\underset{i=1}{\overset{M}{}}\beta _{m(i)}H_{m(i)}(q^{[i]},p^{[i]})\right\},\hfill \\ & =\mathrm{exp}\left\{\underset{m=1}{\overset{M}{}}\beta _mH_m(q^{[i(m)]},p^{[i(m)]})\right\},\hfill \end{array}$$
(16)
where $`i(m)`$ and $`m(i)`$ are the permutation functions in Eq. (3). Then the same derivation that led to the original replica-exchange criterion follows, and the transition probability of replica exchange is given by Eq. (13), where we now have (see Eq. (11))

$$\mathrm{\Delta }=\beta _m\left(E_{\lambda _m}\left(q^{[j]}\right)E_{\lambda _m}\left(q^{[i]}\right)\right)\beta _n\left(E_{\lambda _n}\left(q^{[j]}\right)E_{\lambda _n}\left(q^{[i]}\right)\right).$$
(17)
Here, $`E_{\lambda _m}`$ and $`E_{\lambda _n}`$ are the total potential energies (see Eq. (15)). Note that we need to newly evaluate the potential energy for exchanged coordinates, $`E_{\lambda _m}(q^{[j]})`$ and $`E_{\lambda _n}(q^{[i]})`$, because $`E_{\lambda _m}`$ and $`E_{\lambda _n}`$ are in general different functions.

For obtaining the canonical distributions, the weighted histogram analysis method (WHAM)  is particularly suitable. Suppose we have made a single run of the present replica-exchange simulation with $`M`$ replicas that correspond to $`M`$ different parameter sets $`\mathrm{\Lambda }_m(T_m,\lambda _m)`$ ($`m=1,\mathrm{},M`$). Let $`N_m(E_0,V)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained for the $`m`$-th parameter set $`\mathrm{\Lambda }_m`$. The WHAM equations that yield the canonical probability distribution $`P_{T,\lambda }(E_0,V)`$ with any potential-energy parameter value $`\lambda `$ at any temperature $`T=1/k_B\beta `$ are then given by

$$P_{T,\lambda }(E_0,V)=\left[\frac{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1N_m(E_0,V)}{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1n_me^{f_m\beta _mE_{\lambda _m}}}\right]e^{\beta E_\lambda },$$
(18)
and

$$e^{f_m}=\underset{E_0,V}{}P_{T_m,\lambda _m}(E_0,V).$$
(19)
Here, $`g_m=1+2\tau _m`$, and $`\tau _m`$ is the integrated autocorrelation time at temperature $`T_m`$ with the parameter value $`\lambda _m`$. Note that the unnormalized probability distribution $`P_{T,\lambda }(E_0,V)`$ and the “dimensionless” Helmholtz free energy $`f_m`$ in Eqs. (18) and (19) are solved self-consistently by iteration .

We can use this new replica-exchange method for free energy calculations. We first describe the free-energy perturbation case. The potential energy is given by

$$E_\lambda (q)=E_I(q)+\lambda \left(E_F(q)E_I(q)\right),$$
(20)
where $`E_I`$ and $`E_F`$ are the potential energy for a “wild-type” molecule and a “mutated” molecule, respectively. Note that this equation has the same form as Eq. (15).

Our replica-exchange simulation is performed for $`M`$ replicas with $`M`$ different values of the parameters $`\mathrm{\Lambda }_m=(T_m,\lambda _m)`$. Since $`E_{\lambda =0}(q)=E_I(q)`$ and $`E_{\lambda =1}(q)=E_F(q)`$, we should choose enough $`\lambda _m`$ values distributed in the range between 0 and 1 so that we may have sufficient replica exchanges. From the simulation, $`M`$ histograms $`N_m(E_I,E_FE_I)`$, or equivalently $`N_m(E_I,E_F)`$, are obtained. The Helmholtz free energy difference of “mutation” at temperature $`T`$, $`\mathrm{\Delta }FF_{\lambda =1}F_{\lambda =0}`$, can then be calculated from

$$\mathrm{exp}(\beta \mathrm{\Delta }F)=\frac{Z_{T,\lambda =1}}{Z_{T,\lambda =0}}=\frac{{\displaystyle \underset{E_I,E_F}{}}P_{T,\lambda =1}(E_I,E_F)}{{\displaystyle \underset{E_I,E_F}{}}P_{T,\lambda =0}(E_I,E_F)},$$
(21)
where $`P_{T,\lambda }(E_I,E_F)`$ are obtained from the WHAM equations of Eqs. (18) and (19).

We now describe another free energy calculations based on MREM applied to umbrella sampling , which we refer to as replica-exchange umbrella sampling (REUS). The potential energy is a generalization of Eq. (15) and is given by

$$E_𝝀(q)=E_0(q)+\underset{\mathrm{}=1}{\overset{L}{}}\lambda ^{(\mathrm{})}V_{\mathrm{}}(q),$$
(22)
where $`E_0(q)`$ is the original unbiased potential, $`V_{\mathrm{}}(q)`$ ($`\mathrm{}=1,\mathrm{},L`$) are the biasing (umbrella) potentials, and $`\lambda ^{(\mathrm{})}`$ are the corresponding coupling constants ($`𝝀=(\lambda ^{(1)},\mathrm{},\lambda ^{(L)})`$). Introducing a “reaction coordinate” $`\xi `$, the umbrella potentials are usually written as harmonic restraints:

$$V_{\mathrm{}}(q)=k_{\mathrm{}}\left[\xi (q)d_{\mathrm{}}\right]^2,(\mathrm{}=1,\mathrm{},L),$$
(23)
where $`d_{\mathrm{}}`$ are the midpoints and $`k_{\mathrm{}}`$ are the strengths of the restraining potentials. We prepare $`M`$ replicas with $`M`$ different values of the parameters $`𝚲_m=(T_m,𝝀_m)`$, and the replica-exchange simulation is performed. Since the umbrella potentials $`V_{\mathrm{}}(q)`$ in Eq. (23) are all functions of the reaction coordinate $`\xi `$ only, we can take the histogram $`N_m(E_0,\xi )`$ instead of $`N_m(E_0,V_1,\mathrm{},V_L)`$. The WHAM equations of Eqs. (18) and (19) can then be written as

$$P_{T,𝝀}(E_0,\xi )=\left[\frac{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1N_m(E_0,\xi )}{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1n_me^{f_m\beta _mE_{𝝀_m}}}\right]e^{\beta E_𝝀},$$
(24)
and

$$e^{f_m}=\underset{E_0,\xi }{}P_{T_m,𝝀_m}(E_0,\xi ).$$
(25)
The expectation value of a physical quantity $`A`$ is now given by

$$<A>_{T,𝝀}=\frac{{\displaystyle \underset{E_0,\xi }{}}A(E_0,\xi )P_{T,𝝀}(E_0,\xi )}{{\displaystyle \underset{E_0,\xi }{}}P_{T,𝝀}(E_0,\xi )}.$$
(26)
The potential of mean force (PMF), or free energy as a function of the reaction coordinate, of the original, unbiased system at temperature $`T`$ is given by

$$𝒲_{T,𝝀=\{0\}}(\xi )=k_BT\mathrm{ln}\left[\underset{E_0}{}P_{T,𝝀=\{0\}}(E_0,\xi )\right],$$
(27)
where $`\{0\}=(0,\mathrm{},0)`$.

### 2.2 Replica-Exchange Multicanonical Algorithm

We first briefly review the multicanonical algorithm . Because the coordinates $`q`$ and momenta $`p`$ are decoupled in Eq. (1), we can suppress the kinetic energy part and can write the Boltzmann factor as

$$W_T(x)=W_T(E)=e^{\beta E}.$$
(28)
The canonical probability distribution of potential energy $`P_T(E)`$ is then given by the product of the density of states $`n(E)`$ and the Boltzmann weight factor $`W_T(E)`$:

$$P_T(E)n(E)W_T(E).$$
(29)
In the multicanonical ensemble (MUCA) , on the other hand, each state is weighted by a non-Boltzmann weight factor $`W_{\mathrm{mu}}(E)`$ (which we refer to as the multicanonical weight factor) so that a uniform energy distribution $`P_{\mathrm{mu}}(E)`$ is obtained:

$$P_{\mathrm{mu}}(E)n(E)W_{\mathrm{mu}}(E)\mathrm{constant}.$$
(30)
The flat distribution implies that a free random walk in the potential energy space is realized in this ensemble. This allows the simulation to escape from any local minimum-energy states and to sample the configurational space much more widely than the conventional canonical MC or MD methods.

From the definition in Eq. (30) the multicanonical weight factor is inversely proportional to the density of states, and we can write it as follows:

$$W_{\mathrm{mu}}(E)e^{\beta _0E_{\mathrm{mu}}(E;T_0)}=\frac{1}{n(E)},$$
(31)
where we have chosen an arbitrary reference temperature, $`T_0=1/k_\mathrm{B}\beta _0`$, and the “multicanonical potential energy” is defined by

$$E_{\mathrm{mu}}(E;T_0)=k_\mathrm{B}T_0\mathrm{ln}n(E)=T_0S(E).$$
(32)
Here, $`S(E)`$ is the entropy in the microcanonical ensemble.

A multicanonical Monte Carlo simulation is performed, for instance, with the usual Metropolis criterion: The transition probability of state $`x`$ with potential energy $`E`$ to state $`x^{}`$ with potential energy $`E^{}`$ is given by

$$w(xx^{})=\{\begin{array}{cc}1,\hfill & \mathrm{for}\mathrm{\Delta }E_{\mathrm{mu}}0,\hfill \\ \mathrm{exp}\left(\beta _0\mathrm{\Delta }E_{\mathrm{mu}}\right),\hfill & \mathrm{for}\mathrm{\Delta }E_{\mathrm{mu}}>0,\hfill \end{array}$$
(33)
where

$$\mathrm{\Delta }E_{\mathrm{mu}}E_{\mathrm{mu}}(E^{};T_0)E_{\mathrm{mu}}(E;T_0).$$
(34)
The molecular dynamics algorithm in multicanonical ensemble also naturally follows from Eq. (31), in which the regular constant temperature molecular dynamics simulation (with $`T=T_0`$) is performed by solving the following modified Newton equation:

$$\dot{𝒑}_k=\frac{E_{\mathrm{mu}}(E;T_0)}{𝒒_k}=\frac{E_{\mathrm{mu}}(E;T_0)}{E}𝒇_k,$$
(35)
where $`𝒇_k`$ is the usual force acting on the $`k`$-th atom ($`k=1,\mathrm{},N`$). From Eq. (32) this equation can be rewritten as

$$\dot{𝒑}_k=\frac{T_0}{T(E)}𝒇_k,$$
(36)
where the following thermodynamic relation gives the definition of the “effective temperature” $`T(E)`$:

$$\frac{S(E)}{E}|_{E=E_a}=\frac{1}{T(E_a)},$$
(37)
with

$$E_a=<E>_{T(E_a)}.$$
(38)
The multicanonical weight factor is usually determined by iterations of short trial simulations. The details of this process are described, for instance, in Refs. . However, the iterative process can be non-trivial and very tedius for complex systems.

After the optimal multicanonical weight factor is determined, one performs a long multicanonical simulation once. By monitoring the potential energy throughout the simulation, one can find the global-minimum-energy state. Moreover, by using the obtained histogram $`N_{\mathrm{mu}}(E)`$ of the potential energy distribution $`P_{\mathrm{mu}}(E)`$, the expectation value of a physical quantity $`A`$ at any temperature $`T=1/k_\mathrm{B}\beta `$ can be calculated from

$$<A>_T=\frac{{\displaystyle \underset{E}{}}A(E)n(E)e^{\beta E}}{{\displaystyle \underset{E}{}}n(E)e^{\beta E}},$$
(39)
where the best estimate of the density of states $`n(E)`$ is given by the single-histogram reweighting techniques (see Eq. (30)) :

$$n(E)=\frac{N_{\mathrm{mu}}(E)}{W_{\mathrm{mu}}(E)}.$$
(40)
The replica-exchange multicanonical algorithm (REMUCA)  overcomes both the difficulties of MUCA (the multicanonical weight factor determination is non-trivial) and REM (a lot of replicas, or computation time, is required). In REMUCA we first perform a short REM simulation (with $`M`$ replicas) to determine the multicanonical weight factor and then perform with this weight factor a regular multicanonical simulation with high statistics. The first step is accomplished by the weighted histogram analysis method . Let $`N_m(E)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained at temperature $`T_m=1/k_\mathrm{B}\beta _m`$ of the REM run. The density of states $`n(E)`$ is then given by solving the following WHAM equations :

$$n(E)=\frac{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1N_m(E)}{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1n_me^{f_m\beta _mE}},$$
(41)
where

$$e^{f_m}=\underset{E}{}n(E)e^{\beta _mE}.$$
(42)
Once the estimate of the density of states is obtained, the multicanonical weight factor can be directly determined from Eq. (31) (see also Eq. (32)). Actually, the multicanonical potential energy $`E_{\mathrm{mu}}(E;T_0)`$ thus determined is only reliable in the following range:

$$E_1EE_M,$$
(43)
where

$$\{\begin{array}{cc}\hfill E_1& =<E>_{T_1},\hfill \\ \hfill E_M& =<E>_{T_M},\hfill \end{array}$$
(44)
and $`T_1`$ and $`T_M`$ are respectively the lowest and the highest temperatures used in the REM run. Outside this range we extrapolate the multicanonical potential energy linearly:

$$_{\mathrm{mu}}^{\{0\}}(E)\{\begin{array}{cc}\frac{E_{\mathrm{mu}}(E;T_0)}{E}|_{E=E_1}(EE_1)+E_{\mathrm{mu}}(E_1;T_0),\mathrm{for}E\hfill & <E_1,\hfill \\ E_{\mathrm{mu}}(E;T_0),\mathrm{for}E_1E\hfill & E_M,\hfill \\ \frac{E_{\mathrm{mu}}(E;T_0)}{E}|_{E=E_M}(EE_M)+E_{\mathrm{mu}}(E_M;T_0),\mathrm{for}E\hfill & >E_M.\hfill \end{array}$$
(45)
A long multicanonical MD run is then performed by solving the Newton equations in Eq. (35) into which we substitute $`_{mu}^{\{0\}}(E)`$ of Eq. (45). Finally, the results are analyzed by the single-histogram reweighting techniques as described in Eq. (40) (and Eq. (39)). We remark that our multicanonical MD simulation here actually results in a canonical simulation at $`T=T_1`$ for $`E<E_1`$, a multicanonical simulation for $`E_1EE_M`$, and a canonical simulation at $`T=T_M`$ for $`E>E_M`$ (a detailed discussion on this point is given in Ref. ). Note also that the above arguments are independent of the value of $`T_0`$, and we will get the same results, regardless of its value.

Since the WHAM equations are based on histograms, the density of states $`n(E)`$, or the multicanonical potential energy $`E_{\mathrm{mu}}(E;T_0)`$, will be given in discrete values of the potential energy $`E`$. For multicanonical MD simulations, however, we need the derivative of $`E_{\mathrm{mu}}(E;T_0)`$ with respect to $`E`$ (see Eq. (35)). We thus introduce some smooth function to fit the data. It is best to fit the derivative $`\frac{E_{\mathrm{mu}}(E;T_0)}{E}`$ directly rather than $`E_{\mathrm{mu}}(E;T_0)`$ itself. For this we recall the Newton equation of Eq. (36) and the thermodynamic relation of Eqs. (37) and (38). The effective temperature $`T(E)`$, or the derivative $`\frac{E_{\mathrm{mu}}(E;T_0)}{E}`$, can be obtained by fitting the inverse of Eq. (38) by a smooth function, where $`<E>_{T(E)}`$ is calculated from Eq. (39) by solving the WHAM equations of Eqs. (41) and (42). Given its derivative, the multicanonical potential energy can be obtained by integration:

$$E_{\mathrm{mu}}(E;T_0)=T_0_{E_1}^E\frac{S(E)}{E}𝑑E=T_0_{E_1}^E\frac{dE}{T(E)}.$$
(46)
We remark that the same equations were used to obtain the multicanonical weight factor in Ref. , where $`<E>_{T(E)}`$ was estimated by simulated annealing instead of REM.

Finally, although we did not find any difficulty in the case of protein systems that we studied, a single REM run in general may not be able to give an accurate estimate of the density of states (like in the case of a strong first-order phase transition ). In such a case we can still greatly simplify the process of the multicanonical weight factor determination by combining the present method with the previous iterative methods .

The formulation of REMUCA is simple and straightforward, but the numerical improvement is great, because the weight factor determination for MUCA becomes very difficult by the usual iterative processes for complex systems.

## 3 RESULTS

We now present some examples of the simulation results by the algorithms described in the previous section. Short peptide systems were considered.

For molecular dynamics simulations, the force-field parameters were taken from the all-atom versions of AMBER . The computer code developed in Refs. , which is based on PRESTO , was used. The unit time step was set to 0.5 fs. The temperature during the canonical MD simulations was controlled by the constraint method . Besides gas phase simulations, we have also performed MD simulations with a distance-dependent dielectric, $`ϵ=r`$, and with explicit water molecules of TIP3P model .

As described in detail in the previous section, in generalized-ensemble simulations and subsequent analyses of the data, potential energy distributions have to be taken as histograms. For the bin size of these histograms, we used the values ranging from 0.5 to 2 kcal/mol, depending on the system studied.

The first example is a penta peptide, Met-enkephalin, whose amino-acid sequence is: Tyr-Gly-Gly-Phe-Met. This peptide in gas phase was studied with the force field of AMBER in Ref.  by the replica-exchange MD simulation . We made an MD simulation of $`2\times 10^6`$ time steps (or, 1.0 ns) for each replica, starting from an extended conformation. We used the following eight temperatures: 700, 585, 489, 409, 342, 286, 239, and 200 K, which are distributed exponentially, following the annealing schedule of simulated annealing simulations . As is shown below, this choice already gave an optimal temperature distribution. The replica exchange was tried every 10 fs, and the data were stored just before the replica exchange for later analyses.

As for expectation values of physical quantities at various temperatures, we used the weighted histogram analysis method of Eqs. (41) and (42). We remark that for biomolecular systems the integrated autocorrelation times $`\tau _m`$ in the reweighting formulae (see Eq. (41)) can safely be set to be a constant , and we do so throughout the analyses in this section.

In Figure 1 the time series of temperature exchange (a) and the total potential energy (b) for one of the replicas are shown. We do observe random walks in both temperature space and potential energy space. Note that there is a strong correlation between the behaviors in Figures 1(a) and 1(b).

In Figure 2 the canonical probability distributions obtained at the chosen eight temperatures from the replica-exchange simulation are shown. We see that there are enough overlaps between all pairs of distributions, indicating that there will be sufficient numbers of replica exchanges between pairs of replicas.

We further compare the results of the replica-exchange simulation with those of a single canonical MD simulation (of 1 ns) at the corresponding temperatures. In Figure 3 we compare the distributions of a pair of main-chain dihedral angles $`(\varphi ,\psi )`$ of Gly-2 at two extreme temperatures ($`T=200`$ K and 700 K). While the results at $`T=200`$ K from the regular canonical simulation are localized with only one dominant peak, those from the replica-exchange simulation have several peaks (compare Figures 3(a) and 3(b)). Hence, the replica-exchange run samples much broader configurational space than the conventional canonical run at low temperatures. The results at $`T=700`$ K (Figures 3(c) and 3(d)), on the other hand, are similar, implying that a regular canonical simulation can give accurate thermodynamic quantities at high temperatures.

In Figure 4 we show the average total potential energy as a function of temperature. As expected from the results of Figure 3, we observe that the canonical simulations at low temperatures got trapped in states of energy local minima, resulting in the discrepancies in average values between the results from the canonical simulations and those from the replica-exchange simulation.

We now present the results of replica-exchange umbrella sampling (REUS) . The system of a blocked peptide, alanine-trimer, was studied. Since the thermodynamic behavior of this peptide was extensively studied by the conventional umbrella sampling , it is a good test case to examine the effectiveness of the new method. The force field parameters were taken from the all-atom version of AMBER  with a distance-dependent dielectric, $`ϵ=r`$, which mimics the presence of solvent. We made an MD simulation of $`4\times 10^6`$ time steps (or, 2.0 ns) for each replica, starting from an extended conformation. The data were stored every 20 steps (or, 10 fs) for a total of $`2\times 10^5`$ snapshots.

In Table 1 we summarize the parameters characterizing the replicas for the simulations performed. They are one original replica-exchange simulation (REM1), two replica-exchange umbrella sampling simulations (REUS1 and REUS2), and two conventional umbrella sampling simulations (US1 and US2).

The purpose of the present simulations is to test the effectiveness of the replica-exchange umbrella sampling with respect to the conventional umbrella sampling (REUS1 and REUS2 versus US1 and US2). The original replica-exchange simulation without umbrella potentials (REM1) was also made to set a reference standard for comparison. For REM1 replica exchange was tried every 20 time steps (or, 10 fs), as in our previous work . For REUS simulations, on the other hand, replica exchange was tried every 400 steps (or, 200 fs), which is less frequent than in REM1. This is because we wanted to ensure sufficient time for system relaxation after $`\lambda `$-parameter exchange.

In REM1 there are 16 replicas with 16 different temperatures listed in Table 1. The temperatures are again distributed exponentially. After every 10 fs of parallel MD simulations, eight pairs of replicas corresponding to neighboring temperatures were simultaneously exchanged, and the pairing was alternated between the two possible choices .

For umbrella potentials, the O1 to H5 hydrogen-bonding distance, or “end-to-end distance,” was chosen as the reaction coordinate $`\xi `$ and the harmonic restraining potentials of $`\xi `$ in Eq. (23) were imposed. The force constants, $`k_{\mathrm{}}`$, and the midpoint positions, $`d_{\mathrm{}}`$, are listed in Table 1.

In REUS1 and US1, 14 replicas were simulated with the same set of umbrella potentials at $`T=300`$ K. Let us order the umbrella potentials, $`V_{\mathrm{}}`$ in Eq. (22), in the increasing order of the midpoint value $`d_{\mathrm{}}`$, i.e., the same order that appears in Table 1. We prepared replicas so that the potential energy for each replica includes exactly one umbrella potential (here, we have $`M=L=14`$). Namely, in Eq. (22) for $`𝝀=𝝀_m`$ we set

$$\lambda _m^{(\mathrm{})}=\delta _{\mathrm{},m},$$
(47)
where $`\delta _{k,l}`$ is Kronecker’s delta function, and we have

$$E_{𝝀_m}(q^{[i]})=E_0(q^{[i]})+V_m(q^{[i]}).$$
(48)
The difference between REUS1 and US1 is whether replica exchange is performed or not during the parallel MD simulations. In REUS1 seven pairs of replicas corresponding to “neighboring” umbrella potentials, $`V_m`$ and $`V_{m+1}`$, were simultaneously exchanged after every 200 fs of parallel MD simulations, and the pairing was alternated between the two possible choices. (Other pairings will have much smaller acceptance ratios of replica exchange.) The acceptance criterion for replica exchange is given by Eq. (13), where Eq. (17) now reads (with the fixed inverse temperature $`\beta =1/300k_B`$)

$$\mathrm{\Delta }=\beta \left(V_m\left(q^{[j]}\right)V_m\left(q^{[i]}\right)V_{m+1}\left(q^{[j]}\right)+V_{m+1}\left(q^{[i]}\right)\right),$$
(49)
where replica $`i`$ and $`j`$ respectively have umbrella potentials $`V_m`$ and $`V_{m+1}`$ before the exchange.

In REUS2 and US2, 16 replicas were simulated at four different temperatures with four different restraining potentials (there are $`L=4`$ umbrella potentials at $`N_T=4`$ temperatures, making the total number of replicas $`M=N_T\times L=16`$; see Table 1). We can introduce the following labeling for the parameters characterizing the replicas:

$$\begin{array}{cc}\hfill 𝚲_m=(T_m,𝝀_m)& 𝚲_{I,J}=(T_I,𝝀_J).\hfill \\ \hfill (m=1,\mathrm{},M)& (I=1,\mathrm{},N_T,J=1,\mathrm{},L)\hfill \end{array}$$
(50)
The potential energy is given by Eq. (48) with the replacement: $`mJ`$. Let us again order the umbrella potentials, $`V_J`$, and the temperatures, $`T_I`$, in the same order that appear in Table 1. The difference between REUS2 and US2 is whether replica exchange is performed or not during the MD simulations. In REUS2 we performed the following replica-exchange processes alternately after every 200 fs of parallel MD simulations:

1. Exchange pairs of replicas corresponding to neighboring temperatures, $`T_I`$ and $`T_{I+1}`$ (i.e., exchange replicas $`i`$ and $`j`$ that respectively correspond to parameters $`𝚲_{I,J}`$ and $`𝚲_{I+1,J}`$). (We refer to this process as $`T`$-exchange.)
2. Exchange pairs of replicas corresponding to “neighboring” umbrella potentials, $`V_J`$ and $`V_{J+1}`$ (i.e., exchange replicas $`i`$ and $`j`$ that respectively correspond to parameters $`𝚲_{I,J}`$ and $`𝚲_{I,J+1}`$). (We refer to this process as $`\lambda `$-exchange.)

In each of the above processes, two pairs of replicas were simultaneously exchanged, and the pairing was further alternated between the two possibilities. The acceptance criterion for these replica exchanges is given by Eq. (13), where Eq. (17) now reads

$$\mathrm{\Delta }=\left(\beta _I\beta _{I+1}\right)\left(E_0\left(q^{[j]}\right)+V_J\left(q^{[j]}\right)E_0\left(q^{[i]}\right)V_J\left(q^{[i]}\right)\right),$$
(51)
for $`T`$-exchange, and

$$\mathrm{\Delta }=\beta _I\left(V_J\left(q^{[j]}\right)V_J\left(q^{[i]}\right)V_{J+1}\left(q^{[j]}\right)+V_{J+1}\left(q^{[i]}\right)\right),$$
(52)
for $`\lambda `$-exchange. By this procedure, the random walk in the reaction coordinate space as well as in temperature space can be realized.

We now give the details of the results obtained. In order to confirm that our REM simulations performed properly, we have to examine the time series of various quantites and observe random walks. For instance, the time series of temperature exchange for one of the replicas is shown in Figure 5(a). The corresponding time series of the reaction coordinate $`\xi `$, the distance between atoms O1 and H5, for the same replica is shown in Figure 5(b).

We see that the conformational sampling along the reaction coordinate is significantly enhanced. In the blocked alanine-trimer, the reaction coordinate $`\xi `$ can be classified into three regions : the helical region ($`\xi <3`$ Å), the turn region (3 Å $`<\xi <`$ 7 Å), and the extended region ($`\xi >`$ 7 Å). Thus, Figure 5(b) implies that helix-coil transitions frequently occurred during the replica-exchange simulation, whereas in the conventional canonical simulations such a frequent folding and unfolding process cannot be seen.

After confirming that the present REUS simulations performed properly, we now present and compare the physical quantities calculated by these simulations. In Figure 6 the potentials of mean force (PMF) of the unbiased system along the reaction coordinate $`\xi `$ at $`T=300`$ K are shown. The results are from REM1, REUS1, and US1 simulations. For these calculations, the WHAM equations of Eqs. (24) and (25) were solved by iteration first, and then Eq. (27) was used to obtain the PMF.

From Figure 6 we see that the PMF curves obtained by REM1 and REUS1 are essentially identical for low values of $`\xi `$ ($`\xi <7`$ Å). The two PMF curves start deviating slightly, as $`\xi `$ gets larger, and for $`\xi >9`$ Å the agreement completely deteriorates. The disagreement comes from the facts that the average $`\xi `$ at the highest temperature in REM1 ($`T_{16}=1500`$ K) is $`<\xi >_{T_{16}}8.0`$ Å and that the original REM with $`T`$-exchange only cannot sample accurately the region where $`\xi `$ is much larger than $`<\xi >_{T_{16}}`$. These two simulations were performed under very different conditions: One was run at different temperatures without restraining potentials and the other at one temperature with many restraining potentials (see Table 1). We thus consider the results to be quite reliable for ($`\xi <9`$ Å).

On the other hand, the PMF obtained by US1 is relatively larger than those obtained by REM1 and REUS1 in the region of 2 Å $`<\xi <`$ 4 Å, which corresponds to the structural transition state between the $`\alpha `$-helical and turn structures. This suggests that US1 got trapped in states of energy local minima at $`T=300`$ K. In the region of completely extended structures ($`\xi >9`$ Å), the results of REUS1 and US1 are similar but the discrepancy is again non-negligible. We remark that at $`T=300`$ K the PMF is the lowest for $`\xi 2`$ Å, which implies that the $`\alpha `$-helical structure is favored at this temperature.

We next study the temperature dependence of physical quantities obtained from the REM1, REUS2, and US2 simulations. In Figure 7(a) we show the PMF again at $`T=300`$ K. We observe that the PMF curves from REM1 and REUS2 are essentially identical for $`\xi <9`$ Å and that they deviate for $`\xi >9`$ Å, because the results for REM1 is not reliable in this region as noted above. In fact, by comparing Figures 6 and 7(a), we find that the PMF obtained from REUS1 and REUS2 are almost in complete agreement at $`T=300`$ K in the entire range of $`\xi `$ values shown. On the other hand, we observe a discrepancy between REUS2 and US2 results. The PMF curve for US2 is significantly less than that for REUS2 in the region 2 Å $`<\xi <`$ 8 Å. Note that the PMF curves for US1 and US2 are completely in disagreement (compare Figures 6 and 7(a)).

In Figure 7(b) we show the PMF at $`T=500`$ K, which we obtained from REM1, REUS2, and US2 simulations. We again observe that the results from REM1 and REUS2 are in good agreement for a wide range of $`\xi `$ values. We find that the results from REM1 do not significantly deteriorate until $`\xi >11`$ Å at $`T=500`$ K, whereas it did start deviating badly for $`\xi >9`$ Å at $`T=300`$ K. The PMF curve for US2 deviates strongly from the REUS2 results for $`\xi >6`$ Å and is much larger than that of REUS2 (and REM1) in this region. We remark that at $`T=500`$ K the PMF is the lowest for $`\xi 6`$ Å and low up to $`\xi 8`$ Å, which implies that extended structures are favored at this temperature.

In Figure 8 we show the average values of the reaction coordinate $`\xi `$ of the unbiased system as a function of temperature. The results are again from the REM1, REUS2, and US2 simulations. The expectation values were calculated from Eq. (26). We find that the average reaction coordinate, or the average end-to-end distance, grows as the temperature is raised, reflecting the unfolding of the peptide upon increased thermal fluctuations. Again we observe an agreement between REM1 and REUS2, whereas the results of US2 deviate.

We now present the results of another example of the multidimensional replica-exchange method. This time we consider NPT ensemble of argon fluids, and exchange not only the temperature but also the pressure values of pairs of replicas during a MC simulation . Namely, suppose we have $`M`$ replicas with $`M`$ different values of temperature and pressure ($`T_m`$,$`P_m`$). The state $`x`$ of replica $`i`$ is characterized by the scaled coordinates $`\stackrel{~}{q}^{[i]}`$ and the volume $`𝒱^{[i]}`$ and its weight is given by

$$W(x)=e^{\beta _m\left(E(\stackrel{~}{q}^{[i]})+P_m𝒱^{[i]}\right)+N\mathrm{ln}𝒱^{[i]}}.$$
(53)
The transition probability of replica exchange is then given by Eq. (13), where we now have

$$\mathrm{\Delta }=\left(\beta _m\beta _n\right)\left(E\left(\stackrel{~}{q}^{[j]}\right)E\left(\stackrel{~}{q}^{[i]}\right)\right)+\left(\beta _mP_m\beta _nP_n\right)\left(𝒱^{[j]}𝒱^{[i]}\right).$$
(54)
We prepared $`M=64`$ replicas with $`N_T=8`$ temperature and $`N_P=8`$ pressure values ($`M=N_T\times N_P`$). We alternately exchanged four pairs of temperatures and four pairs of pressures during the replica-exchange simulation. In Figure 9 the values of the set ($`T_m`$,$`P_m`$) in reduced units are shown as crosses.

The shaded regions in Figure 9 are where the acceptance ratios of replica exchange become low ($`<20`$ %). These regions are those where the replica-exchange method fails due to the existence of first-order phase transitions. The results of Figure 9 suggest that the multidimensional REM enables the simulation to connect regions which cannot be reached by one-dimensional REM simulations with only $`T`$-exchange or $`P`$-exchange.

We now present the results of MD simulations based on replica-exchange multicanonical algorithm (REMUCA) . The Met-enkephalin in gas phase was first studied again. The potential energy is, however, that of AMBER in Ref.  instead of Ref. . In Table 2 we summarize the parameters of the simulations that were performed. As discussed in the previous section, REMUCA consists of two simulations: a short REM simulation (from which the density of states of the system, or the multicanonical weight factor, is determined) and a subsequent production run of MUCA simulation. The former simulation is referred to as REM1 and the latter as MUCA1 in Table 2. Finally, a production run of the original REM simulation was also performed for comparison and it is referred to as REM2 in Table 2.

After the simulation of REM1 is finished, we obtained the density of states $`n(E)`$ by the weighted histogram analysis method of Eqs. (41) and (42). The density of states will give the average values of the potential energy from Eq. (39), and we found

$$\{\begin{array}{cc}\hfill E_1& =<E>_{T_1}=30\mathrm{kcal}/\mathrm{mol},\hfill \\ \hfill E_M& =<E>_{T_M}=195\mathrm{kcal}/\mathrm{mol}.\hfill \end{array}$$
(55)
Then our estimate of the density of states is reliable in the range $`E_1EE_M`$. The multicanonical potential energy $`_{mu}^{\{0\}}(E)`$ was thus determined for the three energy regions ($`E<E_1`$, $`E_1EE_M`$, and $`E>E_M`$) from Eq. (45). Here, we have set the arbitrary reference temperature to be $`T_0=1000`$ K.

After determining the multicanonical weight factor, we carried out a multicanonical MD simulation of $`1\times 10^7`$ steps (or 5 ns) for data collection (MUCA1 in Table 2). In Figure 10 the probability distribution of potential energy obtained by MUCA1 is plotted. It can be seen that a good flat distribution is obtained in the energy region $`E_1EE_M`$. In Figure 10 the canonical probability distributions that were obtained by the reweighting techniques at $`T=T_1=200`$ K and $`T=T_M=1000`$ K are also shown. Comparing these curves with those of MUCA1 in the energy regions $`E<E_1`$ and $`E>E_M`$ in Figure 10, we confirm our claim in the previous section that MUCA1 gives canonical distributions at $`T=T_1`$ for $`E<E_1`$ and at $`T=T_M`$ for $`E>E_M`$, whereas it gives a multicanonical distribution for $`E_1EE_M`$.

In the previous works of multicanonical simulations of Met-enkephalin in gas phase (see, for instance, Refs. ), at least several iterations of trial simulations were required for the multicanonical weight determination. We emphasize that in the present case of REMUCA (REM1), only one simulation was necessary to determine the optimal multicanonical weight factor that can cover the energy region corresponding to temperatures between 200 K and 1000 K.

To check the validity of the canonical-ensemble expectation values calculated by the new algorithms, we compare the average potential energy as a function of temperature in Figure 11. In REM2 we used the multiple-histogram techniques (or WHAM) , whereas the single-histogram method  was used in MUCA1. We can see a perfect coincidence of these quantities between REM2 and MUCA1 in Figure 11.

We have so far presented the results of generalized-ensemble simulations of peptides in gas phase. However, peptides and proteins are usually in aqueous solution. We therefore want to incorporate rigorous solvation effects in our simulations in order to compare with experiments. Met-enkephalin was thus studied by both REM and REMUCA simulations in aqueous solution based on TIP3P water model . The AMBER force field of Ref.  was used. The number of water molecules was 526 and they were placed in a sphere of radius of 16 Å. Thirty-six replicas that correspond to temperatures ranging from 200 K to 700 K were used.

The time series of the total potential energy for one of the replicas is shown in Figure 12. We do observe a random walk in potential energy space, which covers an energy range of as much as 2,500 kcal/mol.

For the REMUCA simulation, the multicanonical potential energy and its derivative were obtained by the weighted histogram analysis method from the results of a short REM simulation (of 100 psec). In Figure 13 the probability distribution obtained by the multicanonical production run of this REMUCA simulatoin is plotted. It can be seen that a good flat distribution is obtained in the wide energy range.

Finally, in Figure 14 we compare the distributions of a pair of main-chain dihedral angles $`(\varphi ,\psi )`$ of Gly-2 and Phe-4 around $`T=300`$ K between gas-phase and in-solution results. While the results in gas phase are well localized and sharp, those in aqueous solution are distributed more broadly. This suggests that the energy landscape in gas phase is much more rugged than in aqueous solution; water considerably smoothes out the landscape. We remark that a similar observation was made earlier in Ref. .

## 4 CONCLUSIONS

In this article we reviewed uses of generalized-ensemble algorithms for free-energy calculations in protein folding.

We introduced two new generalized-ensemble algorithms which are generalizations of the replica-exchange method (REM) (we remark that REM is also referred to as parallel tempering). The first one is the multidimensional replica-exchange method (MREM), with which we showed that the replica-exchange method is not limited to tempering (or temperature exchange) and that we can also exchange parameters in the potential energy. One particular realization of this method is replica-exchange umbrella sampling (REUS) where we perform tempering and/or the exchange of parameters that characterize the umbrella potential. The second method is the replica-exchange multicanonical algorithm (REMUCA), in which we combine the merits of REM and multicanonical algorithm (MUCA).

With these new methods available, we believe that we now have working simulation algorithms which we can use for free-energy calculations in protein folding.

Acknowledgements:
We are grateful to Dr. A. Kitao of Kyoto University for useful discussions and collaboration. Our simulations were performed on the Hitachi and other computers at the Research Center for Computational Science, Okazaki National Research Institutes. This work is supported, in part, by a grant from the Research for the Future Program of the Japan Society for the Promotion of Science (JSPS-RFTF98P01101).