統計数理（2006）
第 54 巻 第 2 号 211–222
c 2006 統計数理研究所


特集「予測と発見」［研究詳解］

拡張アンサンブル法による膜タンパク質の立体構造予測1

小久保 裕功 ・岡本 祐幸

2（受付 2006 年 3 月 23 日；改訂 2006 年 7 月 10 日）

要旨

統計力学的アンサンブルを再現する計算機シミュレーション法にマルコフ連鎖モンテカルロ
法があるが，この手法は普通ボルツマン因子に基づいて状態を発生させ，温度一定のカノニカ
ルアンサンブルを実現する．しかし，系の自由度の数が大きいとき，エネルギー極小状態が無数に存在するとともに，それらの間に高いエネルギー障壁が存在するため，シミュレーション
がそれらエネルギー極小状態に留まってしまって，低温における正確なカノニカルアンサンブ
ルを再現するのは困難である．この困難を克服するために開発されたのが，拡張アンサンブル法と総称されるシミュレーション手法である．拡張アンサンブル法は非ボルツマン因子に基づいて，エネルギー空間上の 1 次元酔歩を実現することによって，エネルギー極小状態に留まるのを避ける．ここでは，拡張アンサンブル法の一つである，レプリカ交換法とその一般化版について解説する．そして，具体例として，レプリカ交換モンテカルロ法を膜タンパク質の立体構造予測問題に適用した結果を紹介する．
キーワード： 拡張アンサンブル法，レプリカ交換法，膜タンパク質，膜貫通へリック
ス，タンパク質の立体構造予測．
1.

はじめに

ある多変数関数の最小値を求めることは多くの分野で必要とされることである．その解を予測するために広く使われる手法としては，いくつかの変数値における関数値を求める（観測する）ことによる，または，解が既知の類似関数を求めることによる，帰納的なアプローチがある．一方，計算機シミュレーションなどによって，演繹的に解を求めようとするのが，本稿のアプローチである．
多自由度の複雑系では，エネルギー極小状態が無数に存在し，それらが高いエネルギー障壁
で隔てられているため，従来の温度一定のカノニカルアンサンブル上のマルコフ連鎖モンテカ
ルロ法（MC）などによる計算機シミュレーションでは，それらのエネルギー極小状態に留まってしまって，定量的で正確な予測を行うのが絶望的に難しくなってしまうという難点がある．
このような計算機シミュレーションの困難を克服する有効な方法として，拡張アンサンブル法（generalized-ensemble algorithm）と総称される手法が存在する（筆者による解説としては，例
えば，Mitsutake et al., 2001; 岡本, 2002 などを参照）．この方法は，非ボルツマン因子による
1
2

Department of Chemistry, University of Houston, Houston, TX 77204, U.S.A.
名古屋大学 大学院理学研究科：〒464–8602 愛知県名古屋市千種区不老町

212

統計数理 第 54 巻

第 2 号 2006

人工の統計集団に基づいており，ポテンシャルエネルギー空間上の一次元酔歩を実現することで，エネルギー極小状態に留まるのを避ける．そして，従来の手法と比べて格段に幅広い状態
空間を探索することを可能とするのである．
本稿では，特に広く使われている拡張アンサンブル法の一つである，レプリカ交換法（replica（モンテカルロ法版については，Hukushima and Nemoto, 1996; 分子動力学
exchange method）
法版については，Sugita and Okamoto, 1999; この手法の歴史及び文献の詳細については，Iba,
2001; 伊庭 他，2005 を参照）とその改良版について解説する．そして，具体例として，レプリカ
交換モンテカルロ法の膜タンパク質の立体構造予測シミュレーションへの適用結果を紹介する．
2.

シミュレーション手法

ここでは，レプリカ交換法の詳しい説明を述べる．
2.1 レプリカ交換法
生体系を表す系として，質量が mk（k = 1, . . . , N ）の N 個の原子の系を考え，それらの座標と運動量がそれぞれ q ≡ {q 1 , . . . , q N } と p ≡ {p1 , . . . , pN } で表されるとする．系のハミルトニアン
H(q, p) は運動エネルギー K(p) とポテンシャルエネルギー E(q) の和で次のように与えられる．
(2.1)

H(q, p) = K(p) + E(q) .



ここで，運動エネルギーは次式のとおりである．
N

(2.2)

K(p) =

pk 2
.
2mk
k=1

温度一定のカノニカルアンサンブルでは，各状態 x ≡ (q, p) は次のボルツマン因子の確率重みで実現される．
(2.3)

WB (x; T ) = e−βH(q,p) .

ここで，β はボルツマン定数 kB と温度 T の積の逆数である（= 1/kB T ）．このとき，温度 T における平均運動エネルギーは次で与えられる．

 
N

(2.4)

K(p)T =

k=1

pk 2
2mk

T

3
= N kB T .
2式（2.1）では座標 q と運動量 p が分離されているので，運動エネルギーについては分配関数
中でガウス積分を実行できて，考慮からはずすことができる．よって，ボルツマン因子は次のように書くことができる．
(2.5)

WB (x; T ) = WB (E; T ) = e−βE .

すると，温度一定のカノニカルアンサンブルでは，ポテンシャルエネルギー E の分布 PB (E; T )は状態密度 n(E) とボルツマン因子 WB (E; T ) の積で以下のように与えられる．
(2.6)

PB (E; T ) ∝ n(E)WB (E; T ) .

状態密度はエネルギーと共に急速に増加する関数であり，ボルツマン因子は指数関数的に減少する関数であるから，カノニカル分布は一般にベル型をしている．
この温度 T におけるカノニカル分布を再現する計算機シミュレーション法としては，メトロ
ポリスのモンテカルロ法が広く使われている．この手法では，マルコフ過程に基づきながら，

213

拡張アンサンブル法による膜タンパク質の立体構造予測

以下のメトロポリス判定に従って状態を発生させて行く．すなわち，ポテンシャルエネルギー
E をもつ状態 x から E  をもつ状態 x への遷移確率は次で与えられる．
(2.7)



w(x → x ) = min 1,

WB (x )WB (x)

= min(1, exp (−β∆E)) .

ここで，∆E は次で定義される．
(2.8)∆E ≡ E(x ) − E(x) = E  − E .

実際の計算では，状態 x から次の状態の候補 x を発生させ，∆E を計算する．もし，∆E ≤ 0（ポテンシャルエネルギーが下がる）なら，この候補 x を採択する．逆に，∆E > 0（ポテンシャ
ルエネルギーが上がる）なら，0 から 1 の間の値をとる一様乱数 r を発生させ，r ≤ exp(−β∆E)なら，候補 x を採択するという手順を踏めばよいことになる．
このメトロポリス法は温度一定のカノニカル分布を再現する手法であるが，生体系のように自由度の数が大きく，エネルギー極小状態が無数に存在する系の場合，熱的揺らぎが小さい低温ではシミュレーションがそれらに留まってしまうため，平衡に達するまで膨大な計算時間を必要とする．
この困難を克服する方法として拡張アンサンブル法が考えられたが，まず，レプリカ交換法を説明する．
レプリカ交換法における系は，M 個の違う温度 Tm（m = 1, . . . , M ）をもつ，互いに相互作用しない M 個の独立なコピー（レプリカ）から成る．それぞれのレプリカには常に一つの温度が対応するとする．よって，レプリカと温度には 1 対 1 対応が成り立ち，レプリカを表すラベルを
i（i = 1, . . . , M ），温度を表すラベルを m（m = 1, . . . , M ）とするとき，i と m は置換関係にある．

i = i(m) ≡ f (m) ,m= m(i) ≡ f −1 (i) .

(2.9)

ここで，f (m) は m の置換であり，f −1 (i) はその逆置換である．
[i(M )][1][M ]この系における「状態」を X = {x[i(1)], . . . , xM
} = {xm(1) , . . . , xm(M ) } で表すことにしよう．
1
すなわち，系の状態 X は各レプリカの状態[i]x[i]m ≡ (q )m(2.10)[i]

により指定される．そして，温度 Tm にあるレプリカ i の状態 xm は N 個の原子の座標 q [i] により決まる．
レプリカ同士は相互作用しないので，状態 X の重み因子はボルツマン因子の積で以下のよ
うに与えられる．

M(2.11)WREM (X) = exp −

i=1

 

βm(i) E(q [i] ) = exp −Mβm E(q [i(m)] ).

m=1

ここで，i(m) と m(i) は式（2.9）の置換関数である．
温度がそれぞれ Tm と Tn に対応するレプリカ i と j を交換することを考える．
(2.12)

[j]

[j]
[i]
X = {. . . , x[i]m , . . . , xn , . . . } −→ X = {. . . , xm , . . . , xn , . . . } .

ここで，i，j ，m，n は式（2.9）の置換関数による関係を持つ．よって，レプリカの交換は以下のような新しい置換関数 f  を導入することになる．



214

統計数理 第 54 巻

(2.13)



第 2 号 2006

i = f (m) −→ j = f  (m) ,j = f (n) −→ i = f  (n) .

このレプリカの交換は次のように書くことができる．
(2.14)

[i]

[j]

[j]

[i]

xm ≡ (q [i] )m −→ xm ≡ (q [j] )m ,xn ≡ (q [j] )n −→ xn ≡ (q [i] )n .

ちなみに，この操作はレプリカがそれぞれ i と j に対応する温度 Tm と Tn を交換することと同等である．よって，式（2.14）は次のように書いても良い．

(2.15)

[i]

[i]

[j]

[j]

xm ≡ (q [i] )m −→ xn ≡ (q [i] )n ,xn ≡ (q [j] )n −→ xm ≡ (q [j] )m .

さて，このレプリカ対の交換の操作は，遷移確率 w(X → X  ) に次の詳細釣り合いの条件を課すことによって，平衡状態への収束が保証される．
(2.16)WREM (X)WREM (X  )w(X → X  ) =w(X  → X) ,Z
Z

ここで，Z は系全体の分配関数である．よって，式（2.11）と式（2.16）から以下が得られる．
(2.17)WREM (X  )= exp{−βm E(q [j] ) − βn E(q [i] ) + βm E(q [i] ) + βn E(q [j] )}WREM (X)= exp{−βm [E(q [j] ) − E(q [i] )] − βn [E(q [i] ) − E(q [j] )]}= exp (−∆) .

ここで，∆ は次で与えられる．
∆ = βm (E(q [j] ) − E(q [i] )) − βn (E(q [j] ) − E(q [i] ))(2.18)= (βm − βn )(E(q [j] ) − E(q [i] )) .

(2.19)



ここで，ラベル i，j ，m，n は交換前の置換関数（式（2.9））で次のように関係している．
(2.20)

i = f (m) ,j = f (n) .



よって，レプリカ対の交換確率は，次のメトロポリスの判定条件で与えられる．
(2.21)

[j]w(X → X ) ≡ w(x[i]m | xn ) = min

WREM (X  )1,
WREM (X)= min(1, exp (−∆)) .

[j]
ここで，2 つ目の表現 w(x[i]（と温度対）をあらわに書いた．
m |xn ) では，交換されるレプリカ対
レプリカ交換シミュレーションは，次の 2 つのステップを交互に繰り返すことで実行される（T1 < T2 < · · · < TM とする）．

（1）各レプリカ i（i = 1, . . . , M ）において，一定温度 Tm（m = 1, . . . , M ）の独立なカノニカルア
ンサンブル上の従来の MC シミュレーションを同時にあるステップ数実行する．
（2）隣接し た温 度（Tm と Tm+1 ）に対 応す るレ プリカ 対（i と j ）を 式（2.21）の遷 移確率[i][j]w(xm | xm+1 ) に従って交換する．

215

拡張アンサンブル法による膜タンパク質の立体構造予測

ステップ 2 では隣り合う温度値に対応するレプリカ対を交換することに注意されたい．これは，レプリカ交換の受け入れ確率が温度差に従って指数関数的に減少するからである（式（2.19）と式（2.21）を参照）．また，ステップ 2 において，レプリカ交換が採択されるたびに，式（2.9）の置換関数が更新されることになる．
レプリカ交換シミュレーションは並列計算機に特に適している．各レプリカを各計算機（CPU）に対応させることによって，計算機間の情報通信量を最小にできる訳である（座標の代わりに温度を交換するだけなので）．これは，シミュレーションの間に，MC ステップ t の関数として置換関数 m(i; t) = f −1 (i; t) を記録していくことになる．
長いレプリカ交換シミュレーションの後，物理量 A の任意の温度 T = 1/kB β におけるカノニカルアンサンブル平均は，次式で与えられる．
(2.22)< A >T =E A(E)PB (E; T )E PB (E; T )=−βE
E A(E)n(E) e.
−βE
E n(E)e

ここで，状態密度 n(E) の最適解は次の 2 つの方程式を逐次的に解くことによって得られる（多ヒストグラム再重法（multiple-histogram reweighting techniques）; Ferrenberg and Swendsen,1989）．
(2.23)

M
m=1 Nm (E)M
fm −βm E
m=1 nm en(E) =
e−fm =(2.24)

,n(E) e−βm E .

E

ここで，式中のヒストグラム Nm (E) は，レプリカごとではなく，温度 Tm = 1/kB βm（m =1, . . . , M ）ごとに得られたポテンシャルエネルギー分布のヒストグラムであり，nm はその全サンプル数である．また，fm は無次元化されたヘルムホルツ自由エネルギーである．上の方程式は例えば次のように逐次的に解かれる．先ず，全ての fm（m = 1, . . . , M ）を 0 にセットする．
次に式（2.23）を使って，n(E) を求める．得られた n(E) を式（2.24）に代入して，新しい fm の値を得る．この操作を n(E) と fm の値が収束するまで繰り返す．
更に，ポテンシャルエネルギーの関数でない場合を含む一般の物理量 A のカノニカルアン
サンブル平均は，次式で与えられる．すなわち，まず，上の式（2.23）と（2.24）を解いて，fm（m = 1, . . . , M ）を求める．そして，それを次の式に代入することによって，平均値が得られる（Mitsutake et al., 2003）．
M
m=1(2.25)< A >T =xm A(xm )

M
m=1

xm

exp(−βE(xm ))Mn=1  exp(f − β E(xm ))exp(−βE(xm ))M
=1 n exp(f − β E(xm )).

ここで，xm についての和は温度 Tm に対応するシミュレーションの途中の（トラジェクトリー
の）状態についての和を意味する．
2.2 多次元レプリカ交換法
ここで，レプリカ交換法の多次元（多変数）への拡張版である，多次元レプリカ交換法（multidimensional replica-exchange method; Sugita et al., 2000）を紹介しよう．この手法の開発に至った重要なアイデアは，従来のレプリカ交換法において，M 個のレプリカが相互作用しない
限り，各レプリカのハミルトニアン H(q, p) は同一である必要がないということである．例えば，ハミルトニアンはあるパラメターに依存し，それぞれのレプリカが違うパラメター値をとって

216

統計数理 第 54 巻

第 2 号 2006

いても良いことになる．このとき，温度 Tm に対応する i 番目のレプリカのハミルトニアンは
次のように書くことができる．
Hm (q [i] , p[i] ) = K(p[i] ) + Eλm (q [i] ) .

(2.26)

すなわち，ポテンシャルエネルギー Eλm はパラメター λm に依存しているとする．
従来のレプリカ交換法では，レプリカ i と温度 Tm が 1 対 1 対応しているのに対し，ここでは，レプリカ i と「パラメターセット」Λm ≡ (Tm , λm ) が 1 対 1 対応している．従来のレプリ
カ交換法のパラメター空間が 1 次元（Λm = Tm ）であるので，この手法はその多次元への拡張版ということができる．これによって，温度空間上ばかりでなく，様々なパラメター空間上の酔歩が実行されて，より効率的な状態空間のサンプリングが実現されると期待される．さて，レプリカ同士は相互作用しないので，状態 X の重み因子は各レプリカに対応するボルツマン因子の積で以下のように与えられる（式（2.11）を参照）．

M(2.27)WMREM (X) = exp −[i]βm(i) Eλm (i) (q ) M= exp −

i=1

βm Eλm (q[i(m)]

)



.

m=1

ここで，i(m) と m(i) は式（2.9）の置換関数である．よって，従来のレプリカ交換の採択判定式を導いたのと同じことをすることによって，レプリカ交換法の遷移確率は式（2.21）で与えられることが分かる．ただし，ここでは ∆ は以下のようになる（式（2.18）参照）．
(2.28)∆ = βm (Eλm (q [j] ) − Eλm (q [i] )) − βn (Eλn (q [j] ) − Eλn (q [i] )) .

一般に，全ポテンシャルエネルギー Eλm と Eλn は違う関数なので，新たに，交換後の Eλm (q [j] )と Eλn (q [i] ) を評価する必要があることに注意されたい．
3.

レプリカ交換法による膜タンパク質の立体構造予測

ここでは，レプリカ交換モンテカルロ法の適用例として，膜タンパク質の立体構造予測の結果を紹介する（Kokubo and Okamoto, 2004a, b, c, d）．
具体的には，（SOSUI, Hirokawa et al., 1998 などの WWW サーバーを使って）膜貫通へリック
スの部分だけを取り出し，その膜貫通へリックスだけのレプリカ交換モンテカルロシミュレー
ションを実行した．すなわち，膜貫通へリックスの主鎖の部分は固定して，へリックスの平行
移動と回転及び側鎖の自由度についてのモンテカルロシミュレーションを行った．ポテンシャ
ルエネルギー関数（力場）としては，CHARMM param19 のものを使った．膜及びその外側の水の存在やタンパク質のへリックス以外の部分の存在を考慮するために，以下のような拘束エネ
ルギーも課した（人工的に取り入れたバネの力により，膜から外へ出たり，へリックス同士が不自然に遠くまで行ってしまわないようにした）．




NH −1

(3.1)

Econstr =k1 θ(ri,i+1 − di,i+1 )[ri,i+1 − di,i+1 ]2

i=1

NH

L
L
L 2
{k2 θ(|ziL − z0L | − dLi )[|zi − z0 | − di ]

+

i=1

+


Cα

U
U
U 2
+ k2 θ(|ziU − z0U | − dUi )[|zi − z0 | − di ] }k3 θ (rCα − dCα ) [rCα − dCα ]2 .

拡張アンサンブル法による膜タンパク質の立体構造予測217

ここで，NH は膜貫通へリックスの総数，θ(x) はステップ関数，k1 ，k2 ，k3 はバネ定数，ri,i+1
は i 番目と i + 1 番目のへリックスの隣り合う末端間の距離，di,i+1 はその許容値，ziL (ziU ) はULU
i 番目のへリックスの下端（上端）の z 座標値，z0L (z0U ) はその標準値，dLi (di ) は zi (zi ) のLUz0 (z0 ) からのずれの許容値，rCα は Cα 原子の原点からの距離，dCα は rCα の許容値である．
ここで，許容値とは，対応するずれの値がそれを越えると，バネの力が働いて元に戻そうとす
るようになる「臨界値」のことである．
最初の例は，2 本の膜貫通ヘリックスからなる，glycophorin A である．へリックスは 18 個のアミノ酸からなり，両方とも（一文字表記で）以下のアミノ酸配列を持つ．TLIIFGVMAGVIGTILLI.
レプリカ数は 13 をとり，以下の 13 個の温度を使った．200, 239, 286, 342, 404, 489, 585, 700, 853,
中のパラメターは以下のようにとった．NH = 2，k1 =1041, 1270, 1548, 1888 K. この場合，式（3.1）
2，k3 = 0.05 kcal/（mol Å2 ），di,i+1 = 20 Å，z0L = −13.35 Å，z0U = +13.35 Å，
k2 = 0.5 kcal/（mol Å ）LU
di = di = 1.0 Å，dCα = 50 Å.
まず，レプリカ交換が効率的に行われたことを確かめるために，4 つの時系列を図 1 に示す．
図 1（a）は一つの温度（200 K）が MC sweep とともにどのレプリカに対応したかを示し，図 1（b）∼ 図 1（d）は一つのレプリカ（レプリカ 6）における，それぞれ，温度，ポテンシャルエネル
ギー，自然の構造からの根二乗平均距離（RMSD; 2 つの構造がどれぐらい似ているかの指標でゼロのとき，2 つの構造が完全に一致する）の MC sweep の関数としての変化である．これらの値は全て最小値と最大値の間を酔歩して，レプリカ交換が効率よくなされたことが分かる．
図 2 には，レプリカ交換シミュレーション中に得られた構造のスナップショットを示す．いろいろな構造がサンプルされていることが分かる．
図 3 にはレプリカ交換シミュレーションで得られた立体構造の主成分分析を行った結果を示す．第一主成分軸と第二主成分軸への射影である．高温では，構造の分布がばらばらであるが，低温では，3 つの構造に分類できることが分かった．

図 1.

glycophorin A のレプリカ交換モンテカルロシミュレーションにおける経時変化．
（a）一つの温度（200 K）に対応するレプリカの経時変化，一つのレプリカ（レプリカ 6）
における，（b）温度，（c）ポテンシャルエネルギー，（d）自然の構造からの根二乗平均距
離，の経時変化．

218

統計数理 第 54 巻

第 2 号 2006

図 2.

glycophorin A のレプリカ交換モンテカルロシミュレーションにおけるスナップショッ
ト．
（a）は初期配置．

図 3.

glycophorin A のレプリカ交換シミュレーションで得られた立体構造の主成分分析．
PCA1，PCA2 はそれぞれ第一主成分軸と第二主成分軸に対応．
（a）200 K，（b）342 K，（c）585 K，（d）1888 K.

図 4 では，実験で求められた自然の構造とレプリカ交換シミュレーションで得られた 3 つの構造（図 3）を比較した．クラスター 2 の構造が自然の構造と一致していることが分かる．また，クラスター 1 とクラスター 3 は同一の構造であることも分かる．後者の違いは，2 本のへリックスの番号の付け方の違いに起因するものであり，低温では 2 つの構造しかとらないことが分かる．そして，自然の構造に対応するクラスター 2 が自由エネルギー（及びポテンシャルエネルギー）の最小状態であることが判明した．つまり，この場合，立体構造予測に完全に成功したと言える．

拡張アンサンブル法による膜タンパク質の立体構造予測

図 4.

219

glycophorin A の自然の構造（a）とレプリカ交換シミュレーションで得られた 3 つの構
造（b），
（c），
（d）の比較（図 3 のクラスター 1，2，3 に対応）．（a1），（a2）などの記号の後の番号は一つの構造を違う方向から見た図であることを示す．

最後に，より複雑な系である，bacteriorhodopsin における立体構造予測の結果を紹介する．すなわち，7 本の膜貫通へリックスのレプリカ交換シミュレーションの結果である．この場合，式（3.1）
中のパラメターは以下のようにとった．NH = 7，k1 = k2 = 1.0 kcal/（mol Å2 ）, k3 = 0.05 kcal/（mol Å2 ）, di,i+1 = 20 Å, z0L = 0.0 Å, z0U = 31.5 Å, dLi = dUi = 2.0 Å, dCα = 100 Å.
図 5 にエネルギー極小状態として得られた構造と自然の構造を比較する．7 本のへリックスの配置が一致していることが分かる．

220

統計数理 第 54 巻

図 5.

4.

第 2 号 2006

bacteriorhodopsin の自然の構造（a）とレプリカ交換シミュレーションで得られたエネ
ルギー極小構造（b）の比較．（a1），（a2）などの記号の後の番号は一つの構造を違う方向から見た図であることを示す．自然の構造（a）では，シミュレーションでは取り入れな
かったレチナール分子を紫で表した．

おわりに

本稿では，タンパク質の立体構造予測問題において有効な演繹的手法として，拡張アンサン
ブル法について解説した．特にレプリカ交換法は従来の温度一定のカノニカルシミュレーショ
ンを並行して実行するだけなので，適用がきわめて容易である．実際，レプリカ交換分子動力学法はタンパク質の折り畳み問題などにおいて広く使われている．

参 考 文 献Ferrenberg, A. M. and Swendsen, R. H.（1989）. Optimized Monte Carlo data analysis, PhysicalReview Letters, 63, 1195–1198.
Hirokawa, T., Boon-Chieng, S. and Mitaku, S.（1998）. SOSUI: Classification and secondary structure
prediction system for membrane proteins, Bioinformatics, 14, 378–379.
Hukushima, K. and Nemoto, K.（1996）. Exchange Monte Carlo method and application to spin glass
simulations, Journal of Physical Society of Japan, 65, 1604–1608.
Iba, Y.（2001）. Extended ensemble Monte Carlo, International Journal of Modern Physics C, 12,
623–656.
伊庭幸人，種村正美，大森裕浩，和合 肇，佐藤整尚，高橋明彦（2005）. 『計算統計 II マルコフ連鎖モンテカルロ法とその周辺』，統計科学のフロンティア 12, 岩波書店，東京．
Kokubo, H. and Okamoto, Y.（2004a）. Prediction of transmembrane helix configurations by replicaexchange simulations, Chemical Physics Letters, 383, 397–402.
Kokubo, H. and Okamoto, Y.（2004b）. Prediction of membrane protein structures by replica-exchange拡張アンサンブル法による膜タンパク質の立体構造予測

221

Monte Carlo simulations: Case of two helices, Journal of Chemical Physics, 120, 10837–10847.
Kokubo, H. and Okamoto, Y.（2004c）. Classification and prediction of low-energy membrane protein
helix configurations by replica-exchange Monte Carlo method, Journal of Physical Society ofJapan, 73, 2571–2585.
Kokubo, H. and Okamoto, Y.（2004d）. Self-assembly of transmembrane helices of bacteriorhodopsin
by a replica-exchange Monte Carlo simulation, Chemical Physics Letters, 392, 168–175.
Mitsutake, A., Sugita, Y. and Okamoto, Y.（2001）. Generalized-ensemble algorithms for molecularsimulations of biopolymers, Biopolymers（Peptide Science）, 60, 96–123.
Mitsutake, A., Sugita, Y. and Okamoto, Y.（2003）. Replica-exchange multicanonical and multicanonical replica-exchange Monte Carlo simulations of peptides. I. Formulation and benchmark test,
Journal of Chemical Physics, 118, 6664–6675.

岡本祐幸（2002）. 拡張アンサンブルの方法，『生体系のコンピュータ・シミュレーション』（岡崎 進，
岡本祐幸 編），化学フロンティア No. 8, 25–42, 化学同人，京都．
Sugita, Y. and Okamoto, Y.（1999）. Replica-exchange molecular dynamics method for protein folding,
Chemical Physics Letters, 314, 141–151.
Sugita, Y., Kitao, A. and Okamoto, Y.（2000）. Multidimensional replica-exchange method for freeenergy calculations, Journal of Chemical Physics, 113, 6042–6051.

222

Proceedings of the Institute of Statistical Mathematics Vol. 54, No. 2, 211–222(2006)

Prediction of Membrane Protein Structuresby Generalized-ensemble Algorithms
Hironori Kokubo1 and Yuko Okamoto2
1 Department of Chemistry, University of Houston
2 Department of Physics, Nagoya University

The Markov-chain Monte Carlo method is a computer simulation algorithm that reproduces statistical ensembles. It is usually based on the Boltzmann weight factor andrealizes a ﬁxed-temperature canonical ensemble. However, when the number of degrees of
freedom of the system is large, there exist a huge number of local-minimum-energy statesthat are separated by high-energy barriers. This forces the simulation to get trapped inenergy-local-minimum states and makes it very diﬃcult to reproduce an accurate lowtemperature canonical ensemble. Generalized-ensemble algorithms are generic terms for
those methods that are based on non-Boltzmann weight factors and overcome the abovementioned diﬃculty by realizing a one-dimensional random walk in energy space. We
review one of the generalized-ensemble algorithms, namely, the replica-exchange method,and its extensions. As an example of its application, we present the results of replicaexchange Monte Carlo simulation applied to the prediction of membrane protein structures.

Key words: Generalized-ensemble algorithm, replica-exchange method, membrane protein, transmembrane helix, protein tertiary structure prediction.