複雑系の	  拡張アンサンブルアルゴリズム	
光武 亜代理（慶應大•理工）	

目的とこれまでの研究	
化学物理に基づく全原子の計算機シミュレーションによって蛋白質のフォール
ディングや会合系などの機構について理解する	
フォールディング、会合系のシミュレーション	
  （たくさんの極少値からなる自由エネルギー曲面）	
  • １つの極少状態にトラップされる	
  -­‐>	
  効率の良いサンプリング方法が必要	  – 拡張アンサンブル法をペプチド系へ適用し有効性を調べた	
  – 複合型拡張アンサンブル法の開発を行った	  （近年、多変数版（会合系シミュレーションのため））	
  – 陽に水を取り入れたProteinGの大規模なフォールディングシミュレーション	  – 分子の有限サイズ系への拡張アンサンブル法の適用	  •正確なエネルギー関数が必要（特に、溶媒を正確に取り入れる）	
  – 陽に水分子を加える	  • 近年、専用ボードMD-­‐GRAPE3やGPGPUを用いてシミュレーションを行っている	  – 統計力学に基づく液体論であるRISM	
  理論を用いる	
  • 1D	
  RISM理論を用いた拡張アンサンブルシミュレーション	
  • ３D	
  RISM理論を用いてペプチドと蛋白質の会合系への適用	
  緩和モード解析手法を蛋白質系へ適用	
  – ペプチド系で有効性を調べた	
  – リゾチウムという蛋白質系に適用	
  タンパク質	n20種類のアミノ酸がペプチド結合により

	  つながった1本鎖からなる高分子	  n１つのアミノ酸配列に対して１つの天然 	
  
立体構造	
  
n73,338個の立体構造（2011年5月30日） 	
  (h?p://www.pdb.org/pdb/home/home.do より）	  n2,000,000個以上の遺伝子（2006年）	

２次構造	
  
αヘリックス（ピンク）	
  βシート（黄色）	

タンパク質系のシミュレーション	
Monte	  Carlo	
  法： 重み因子に基づいて確率論的に状態を発生
分子動力学法： ニュートン方程式を解く	
WB(U)=exp(-­‐βU)	
  m!!
ri = fi = "!U (r )
!ri

共有結合力	

非共有結合力	+溶媒との相互作用エネルギー�Ø厳密に水の効果を取り入れる(分子論的取り扱い）	

統計力学に基づくRISM理論	
水をあらわに取り入れる 		  タンパク質のフォールディング	
•
•
•
•
•

天然条件で1つの構造をもつが温度や溶媒を変えると変性	
  条件をもどすと巻き戻る 〜100残基ぐらい	  小タンパク質(100残基程度）のフォールディング速度はμsからms	
  2次構造形成（ヘリックス） 10μｓ	  
最低 →１μsは必要	
  シミュレーションのΔt=1fs	
  2状態フォールディングの２つの脱出速度	

ku = vu exp(# !$Gu0" )k f = v f exp(# !$G 0f " )kf
0"1#G = " ! ln( )ku
シミュレーション可能な時間では,通常の方法だとエネルギー極小状態に留まってしまうので,エネルギー障壁を越えるサンプリング手法が必要	タンパク質系のプログラムパッケージ	
	
  	
  	
  	
  	
   ECEPP(ECEPP)	

AMBER(AMBER(f94,f96,f99,f03)) http://ambermd.org/	
CHARMM(CHARMM22) http://www.charmm.org/	
NAMD (CHARMM) http://www.ks.uiuc.edu/Research/namd/	
GROMACS(GROMOS) 	
http://www.gromacs.org/Documentation/Terminology/Force̲Fields/GROMOS	TINKER(AMBER(P94,P96,OPLS-AA.CHARMM22)	
IMPACT(AMBER,OPLS) http://lutece.rutgers.edu/levygroup/	
PRESTO(AMBER(P89,P94,P96),OPLS)	
PEACH(AMBER,OPLS) http://staff.aist.go.jp/y-komeiji/index̲j.html	
MARBLE	COSMOS90（AMBER)	
専用ハードウェア （約2万原子）	
MDGRAPE3~数十ns (5ns/日 ）	
GPGPU (NAMD, AMBER) ~数十ns (30ns/日 ）	
ANTON (Shaw group)~数百ns (17μs/日512node） 数十ns	拡張アンサンブル法	
n物理の磁性の分野で開発(MC（人工の重み因子）=> MD（人工のポテンシャル）)	

n人工の重み因子（人工のポテンシャル）を用いてシミュレーションを行なう	
カノニカルモンテカルロ	

ボルツマン因子�

WB(E) = exp(! "E)
小タンパク質の巻き戻り	  μs~ms	
  MDの時間ステップ fs	  拡張アンサンブル法	

人工の重み因子 	

WGEN

エネルギー空間をランダムウォークすることにより広い構造空間の探索を実現	
利点 唯一回のシミュレーションの結果から、最小エネルギー状態
	
ばかりでなく、幅広い温度領域における熱力学量を計算できる	→Reweighting	拡張アンサンブル法の利点１	
nエネルギー空間上のランダムウォークをするため極小エネルギー状態にトラップされずに幅広い構造空間をサンプリングできる	ØエネルギーのTime series 	

低温のカノニカルシミュレーション (T=200K)	

Helix-coil transitions of a peptide	
マルチカノニカルシミュレーション�-60-80

E (kcal/mol)

-100
-120
-140
-160
-180
-200
120000 140000 160000 180000 200000 220000 240000 260000 280000 300000
Nsweep極小状態に	

トラップ�

拡張アンサンブル法の利点2	
  • Reweighting法を用いて任意の温度の熱力学量を再現できる	
• Single Histogram Reweighting法	Ø 1つの拡張アンサンブルシミュレーションから再現	
M% A(x )exp(#$E(x ))W (E(x ))<A> "% exp(#$E(x ))W (E(x ))
i

i=1

T

#1
mu

i

M

i

i=1

#1
mu

i

i

• Multiple Histogram Reweighting 法	!

Ø 多数のシミュレーション結果から再現	
Rn( E ) = R"1j# g H (E)
j =1

"1
j

j

# g h exp(" ! E + f )
j =1

j

j

j

exp(" f j ) = $ dEn( E ) exp(" ! j E ) ! # n( E ) exp(" ! j E )E既存の拡張アンサンブル法	
マルチカノニカル法 (Multicanonical Algorithm)	
���B. A. Berg & T. Neuhaus, Phys. Lett. B267, 249 (1991).B. A. Berg & T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992).	

	

	

	

	

	

Wmu(E) !n-1 (E)Emu ( E; T0 ) = E0 ( E; T0 ) + kbT0 ln P( E )Pmu(E) ! n(E)Wmu(E) = const

シミュレーティドテンパリング (Simulated Tempering)!
A.P. Lyubartsev, et al., J. Chem. Phys. 96, 1776 (1992). E. Marinari and G. Parisi, Europhys. Lett. 19, 451 (1992).	

	
W (E;T ) = exp(! " E + a )
ST
m
m
m
	
P (T) = ! dEn(E) exp(" #E + a(T)) = const	 ST	 exp( !am )"# dEn(E) exp( !$m E)
レプリカ交換法 (Replica-Exchange Method)	
����(PARALLEL TEMPERINGおよびMULTIPLE MARKOV CHAIN METHOD）�	

K.Hukushima & K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996). 	

K. Hukushima, H. Takayama & K. Nemoto, Int. J. Mod. Phys. C 7, 337 (1996).!
M

WREM(X) = ! exp$%&" #iE(x )'()ii =1

マルチカノニカル法	
  (MUCA)	
  
B. A. Berg & T. Neuhaus, Phys. Lett. B267, 249 (1991)."
B. A. Berg & T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992)."
Multicanonical
Canonical 50K
Canonical 1000K

Pmu(E) ! n(E)Wmu(E) = const

Wmu(E) !n-1 (E)

30

E

20
10
0
-10
0

エネルギー空間上で平らな確率分布を実現	

エネルギー空間上のランダムウォークを実現		

重みを決めるために何回かの短いシミュレーショ
ンが必要	
From H {n!1} ( E ) = n( E )Wmu{n!1} ( E )Wmu{n} ( E ) = 1/ n( E ) = Wmu{n!1} ( E ) / H {n!1} ( E )50000

100000MC Sweeps

150000200000ＭＣ−＞ＭＤ（カノニカルアンサンブル(NVT)）	
  
f (q, p ) ! exp(" H (q, p ) / k BT )Nose,	
  S.,Prog.	
  of	
  Theo.	
  Phys.	
  Supp.	
  103,	
  1(1991)	

1. Gauss thermostat	dq i p i
= ,
dt
mi

p i2
= consti =1 2 m
i

! =$

pi
g
" k BT0 ) exp("U (q ) / k BT )
i 2m
2
i

g = 3N ! 1

dp i"U (q )=#
# ! pidt"q i 2p i "UNN

K =!

f ( p, q ) # ! ( $% m # "q
i =1

i

p
%
i =1 2 miN

2. Stochastic method 	dq i p i dp i
#U
= ,
="
" ! pi + R i .
dt
mi dt
#q i

R (t ) = 0,
R! (t ) R" (t %) = 2k BT0#$ !" $ (t & t %)p i2
f (q, p ) = exp(! H (q, p ) / k BT0 ), H (q, p ) = "+ U (q)i 2mi

3 Extended system (Nose-Hoover thermostat)	

dq i p i dp i
#U (q )
= ,="" ! pi ,
dt
mi dt
#q i
p i2
d!
1
g = 3N= ($ " gk BT0 ).
dt Q i mi2pf (q, p ) ! exp(" H 0 (q, p ) / k BT0 ), H 0 (q, p ) = # i + U (q )
i 2m
i

2
i

i

N

=$

p i "U% m # "qi =1i

gk BT0i

Multicanonical MD	

U.H.E. Hansmann, Y. Okamoto and F. Eisenmenger, Chem. Phys. Lett., 259, 321 (1996).	

N. Nakajima, H. Nakamura and A. Kidera, J. Phys. Chem. B, 101, 817 (1997)．	

T. Munakata and S. Oyama, Phys. Prev. E, 54, 4394 (1996).	

U mu (U )マルチカノニカルポテンシャル�We want a simulation which realizes random walk in potential energy
space.	

Pmu (U ) ! n(U )Wmu (U )Pca (U ) ! n(U ) exp("U / k BT0 )= n(U ) exp("U mu (U ) / k BT0 ) = const) でシミュレーション�
温度 T0 でポテンシャルエネルギ U mu (U����!U mu (U (q ))!U (q )"!q i
!q i
!U mu (U ( q ))!U mu (U ) !U !U mu (U ) 0
mu
fi = "
="
i
=
if i
!q i
!U
!q i
!U
!U mu (U )We have to get!U. The process is the same to MC	コメント�
1. アンブレラサンプリングとの対応�
See:C. Bartels and M. Karplus, J. Phys. Chem. B, 102, 865(1998).	

U mu (U ) = k BT0 ln n(U ) = U + k BT0 ln Pc (U )From Pc (U ) = n( E ) exp(!U / k BT0 )2. the scale factorの意味�
Y. Sugita & Y.Okamoto, Chem. Phys. Lett. 329, 261 (2000).	

�

A. Mitsutake, Y. Sugita, & Y.Okamoto., Biopolymers (Peptide Science) 60, 96(2001)	

"U mu (U !)T" ln n(U !)"S (U !)
# k BT0= T0= 0"U ! U !=U"U ! U !=U"U ! U !=U T (U )U =< U >T (U )At low energy T(U)<T0	

T0>1
T (U )

Enhance to change 	

conformation	シミュレーティッド・テンパリング法	
	A.P.	
  Lyubartsev,	
  et	
  al.,	
  J.	
  Chem.	
  Phys.	
  96,	
  1776	
  (1992).	
  	
  
E.	
  Marinari	
  and	  G.	
  Parisi,	
  Europhys.	
  Le4.	
  19,	
  451	  (1992).	
  Protein	
  system:	
  U.H.E	
  Hansmann,	
  Y.	
  Okamoto,	
  J.	
  Compt.	
  Chem.	
  18.	
  920	  (1997)	
  温度が動的変数：	
  	
  温度空間上をランダムウォークするような重みを用いる	
  WST(E;T) = exp(! "E + a(T))
PST (T) = ! dEn(E) exp(" #E + a(T)) = const

exp( !am )"# dEn(E) exp( !$m E)

シミュレーション中に重みに基づき温度を変える	
  シミュレーションの振る舞いはレプリカ交換に似ている	
  温度空間上のランダムウォーク −＞ エネルギー空間上のランダムウォーク	  	
  欠点：重み因子を決める必要があり、系が大きくなると難しい	
  利点：１つのレプリカで良い（レプリカ間の通信が必要ない−＞ＧＰＧＰＵなどに適する？）	
  シミュレーティッド・テンパリング法	
• N個の温度を用意する	
  • シミュレーションの間	
  1.	
  温度Ｔｍでの普通のMC	
  or	
  MD	
  シミュレーションを行なう	
  	
  
2.	
  あるステップおきに温度 Tm を隣の温度 Tm ± 1 に変える	

amを決定するのに短いSTシミュレーションを何度か逐次的に行う必要がある。�
シミュレーティッドテンパリング法は一度重み因子が決定すれば、	

パワフルな方法だが、系が複雑になると重み因子をきめるのが難しい。�

レプリカ交換法	
  (REM)	
  (PARALLEL TEMPERINGおよびMULTIPLE MARKOV CHAIN METHODとも呼ばれる） 	

K. Hukushima & K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996)	

�K. Hukushima, H. Takayama & K. Nemoto, Int. J. Mod. Phys. C 7, 337 (1996).	

1. 系	
異なる温度における互いに相互作用しない系の	
コピー（レプリカ）をＭ個用意	
M
	 WREM(X) = ! exp$%" #iE(xi )'( X = {x1,x 2,..., xM}
i =1

&

)

2. Replica-Exchange	
・各レプリカで独立にカノニカルシミュレーションを実行		%'
1 for " # 0" = ! {( E ( X ') # E ( X )}	
 w(x ! x' )= &'exp($") for " > 0(
	
・温度値が隣接した２つのレプリカを詳細つり合いの条件	
に基づいて交換	" = [ ! n # ! m ]( E (q[i ( m )] ) # E (q[ f ( n )] ))この２つの操作を繰り返す 	
3. 各温度のカノニカル分布を再現	
Multiple Histogram Reweighting法 を用いて	
任意の温度のカノニカル分布計算	
�・レプリカはいろいろな温度の状態をサンプルするので、		
低温で極小状態にトラップされない。→ 低温でのサンプリング効率よい	
・重み因子を決定する必要がない。
並列コンピュータに適している		Replica-­‐exchange	
  Method	
  (REM)	
  (also	
  referred	
  to	
  as	
  Parallel	
  tempering	
  and	  Muliple	  Markov	
  chain	
  method)	
  	
  
�

K.	
  Hukushima	
  &	
  K.	
  Nemoto,	
  J.	
  Phys.	
  Soc.	
  Jpn.	
  65,	
  1604	
  (1996),	
  K.	
  Hukushima,	
  H.	
  Takayama	
  &	
  K.	
  Nemoto,	
  Int.	
  J.	
  Mod.	
  Phys.	
  C	
  7,	
  337	  (1996).	
  M

WREM(X) = ! exp$%&" #iE(x )'() X = {x1,x 2,..., xM}ii =1

Product of boltzman factors at some temperatures	Replica-­‐exchange	
  Method	
  (REM)	
  (also	
  referred	
  to	
  as	
  Parallel	
  tempering	
  and	  Muliple	  Markov	
  chain	
  method)	
  	
  
�

K.	
  Hukushima	
  &	
  K.	
  Nemoto,	
  J.	
  Phys.	
  Soc.	
  Jpn.	
  65,	
  1604	
  (1996),	
  K.	
  Hukushima,	
  H.	
  Takayama	
  &	
  K.	
  Nemoto,	
  Int.	
  J.	
  Mod.	
  Phys.	
  C	
  7,	
  337	  (1996).	
  M

WREM(X) = ! exp$%&" #iE(x )'() X = {x1,x 2,..., xM}ii =1

Product of boltzman factors at some temperatures	

↓	
  

↓	
  

↓	
  

Replica-­‐exchange	
  MD	
  method	
  (REM)	
  (also referred to as Parallel tempering and Multiple Markov chain method) 	

� MD: Y.Sugita and Y. Okamoto, Chem. Phys.Lett. 314, 141 (1994)"

xm[ i ] ! (q[ i ] , p[ i ] ) m"�

Tm (m = 1,! , M )X = {x1[ i (1)] ,! , xM[ i ( M )] } = {xm[1](1) ,! , xm[ M( M] ) }M

T1 < T2 < ! < TM

i: label of replcia	

m: label of temperature	

W REM ( X ) = ! exp(!! m ( i ) H (q [ i ] , p [ i ] )), where H = K ( p [ i ] ) + U (q [ i ] )
i =1

[i ]
m

$ x " (q[ i ] , p[ i ] ) m # xm[ j ]! " (q[ j ] , p[ j ]! ) m% [ j][ j][ j][ i ]!
[i ][ i ]!
& xn " (q , p ) n # xn " (q , p ) n

p[ i ]!

Tn [ i ]
"
p ,
Tm

p

[ j ]!

Tm [ j ]
"
p
Tn

!"0
# 1,
Wn ,m = $" = [ ! n # ! m ]( E (q[i ( m )] ) # E (q[ f ( n )] ))&exp(%!), ! > 0
MCと同じ�（運動量のスケーリングの詳細は	

Mori Y. and Okamoto Y., J. Phys. Soc. Japan, 79, 074001 (2010)）�

レプリカ交換シミュレーションの結果	
T	E	

Low energy	

60000	

Low energy	

High energy	

Once extend	

160000	

260000	

レプリカ交換から�

従来のカノニカル法から�目的とこれまでの研究	
化学物理に基づく全原子の計算機シミュレーションによって蛋白質のフォール
ディングや会合系などの機構について理解する	
フォールディング、会合系のシミュレーション	
  （たくさんの極少値からなる自由エネルギー曲面）	
  • １つの極少状態にトラップされる	
  -­‐>	
  効率の良いサンプリング方法が必要	  – 拡張アンサンブル法をペプチド系へ適用し有効性を調べた	
  – 複合型拡張アンサンブル法の開発を行った	  （近年、多変数版（会合系シミュレーションのため））	
  – 陽に水を取り入れたProteinGの大規模なフォールディングシミュレーション	  – 分子の有限サイズ系への拡張アンサンブル法の適用	  •正確なエネルギー関数が必要（特に、溶媒を正確に取り入れる）	
  – 陽に水分子を加える	  • 近年、専用ボードMD-­‐GRAPE3やGPGPUを用いてシミュレーションを行っている	  – 統計力学に基づく液体論であるRISM	
  理論を用いる	
  • 1D	
  RISM理論を用いた拡張アンサンブルシミュレーション	
  • ３D	
  RISM理論を用いてペプチドと蛋白質の会合系への適用	
  緩和モード解析手法を蛋白質系へ適用	
  – ペプチド系で有効性を調べた	
  – リゾチウムという蛋白質系に適用	
  複合型拡張アンサンブル法の開発	
REVIEW: A. Mitsutake, Y. Sugita, and Y.Okamoto., Biopolymers (Peptide Science) 60, 96 (2001); cond-mat/0012021	
連載：拡張アンサンブル法入門（I-­‐V) 2005年-­‐2006年分子シミュレーション研究会会誌

	

レプリカ交換法（REM)	
  

	
  	
  利点：重み因子を決める必要がない（最も簡単な拡張アンサンブル法）	
  	
  	
  欠点：系が大きくなると、必要なレプリカの数が膨大になる	  従来のマルチカノニカル法(MUCA)とシミュレーティッドテンパリング法(ST)	
  	
  	
  利点：たくさんのレプリカを必要としない	
  

	
  	
  欠点：重み因子の決定の困難（根気と熟練が必要）�

1.	
  レプリカ交換マルチカノニカル法	  (REMUCA)	
  
	
  	
  	
  	
  マルチカノニカルレプリカ交換法	
  (MUCAREM)	
  Y.	
  Sugita	  	
  and	
  Y.Okamoto,	
  Chem.	
  Phys.	
  Le4.	
  329,	
  261	  (2000).	
  
	
  	
  	
  	
  	
  	
  	
  REMUCA	  and	
  MUCAREM	  for	  Monte	  Carlo	  version.	
  
	
  	
  	
  	
  	
  	
  	
  A.Mitsutake,	
  Y.Sugita,	
  and	
  Y.Okamoto,	
  J.Chem.Phys,	
  118	  (2003)	
  6664	
  	
  
	
  	
  	
  	
  	
  	
  	
  A.Mitsutake,	
  Y.Sugita,	
  and	
  Y.Okamoto,	
  J.Chem.Phys,	
  118	  (2003)	
  6676	
  

2.交換シミュレーティッドテンパリング法	
  	  (REST)	
  A.	
  Mitsutake	  	
  and	
  Y.Okamoto,	
  Chem.	
  Phys.	
  Le4.,	
  332,	
  131	  (2000).	
  3.シミュレーティッドテンパリングレプリカ交換法（STREM)	
  A.Mitsutake	
  and	  Y.	
  Okamoto,	
  J.	
  Chem.	
  Phys.	
  121,2491(2004)	
  

4．REST	
  and	
  STREM	
  の多変数版	
  A.	
  Mitsutake	
  and	  Y.	
  Okamoto,	
  J.	
  Chem.	
  Phys.	
  130	  (2009)	
  214105	
  A.	
  Mistsutake	
  and	  Y.	
  Okamoto,	
  Phys.	
  Rev.	
  E	
  79	
  (2009)	
  214105	
  A.	
  Mitsutake,	
  J.	
  Chem.	
  Phys.	
  131	  (2009)	
  094105.	
  	

WMU ( E;Tm ) ! 1/ n( E )M

Wmucarem(E) = ! Wmui(Ei )i =1

WST(E;Tm ) = exp(! "m E + am )M

Wstrem( E ) = !Wst i ( Ei )i=1

短いレプリカ交換シミュレーションから	  重みを決定	

Multiple-histogram reweighting techniques	

A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 63, 1195 (1989).	

n複数のシミュレーション結果を用いて熱力学量を計算する	
Rn( E ) = R"1j# g H (E)
j =1

"1
j

j

# g h exp(" ! E + f )
j =1

j

j

j

exp(" f j ) = # dEn( E ) exp(" ! j E )! $ n( E ) exp(" ! j E )Ehi(E):�温度Tiのエネルギーヒストグラム	

Hi:�温度Tiの全サンプル数	

n(E):�状態数	

fi: 温度Ti�のDimensionless Helmholtz free energy 		

これらの方程式の解は逐次計算により得られる�
n 1/n(E)：マルチカノニカルパラメーターに対応	
  n ｆ（T):シミュレーティッド・パラメーターに対応	  短いレプリカ交換シミュレーションの結果から	
  各アンサンブルシミュレーションのパラメータを求めることができる	

複合型拡張アンサンブル法の結果	
n簡単に重み因子を決めることができた	
->より大きな系へ	
nレプリカ交換法よりもレプリカの数が少ない	n幅広いエネルギー（温度）領域のときは	
MUCAREMまたはSTREMからREMUCAまたはREST	
->より大きな系へ	

REM	

REMUCA	

MUCAREM	-2-3-4

いくつかの	

マルチカノニカル	

アンサンブル	-5ln P(E)

マルチカノニカル	

アンサンブル	

-6
-7
-8
-9
-10

REST	

シミュレーティッド	テンパリング	アンサンブル	-11-220 -200 -180 -160 -140 -120 -100 -80
E (kcal/mol)

-60

-40

STREM	

いくつかの	

シミュレーティッド	テンパリング	アンサンブル	多変数RESTとSTREMの開発	
•A.	
  Mitsutake	
  and	  Y.	
  Okamoto,	
  J.	
  Chem.	
  Phys.	
  130	  (2009)	
  214105	
  A.	
  Mistsutake	
  and	  Y.	
  Okamoto,	
  Phys.	
  Rev.	
  E	
  79	
  (2009)	
  214105	
  A.	
  Mitsutake,	
  J.	
  Chem.	
  Phys.	
  131	  (2009)	
  094105.	
  会合系のシミュレーション	  – バイアスポテンシャル上のランダムウォークも実現	L

エネルギー	E " (x) = E 0 (x) + # "(l ) Vl (x)分配関数（ T	
  and	  λ）	<-­‐-­‐	  d	
  -­‐-­‐>	  l=1

Z(T,") = % exp(#$ E " (x))
xEx.1：Umbrella	
  sampling	  E 0 + "V $'"J # Si S j + h&"# Si )% i (<i, j>! system	  Ex.2：Spin	  !

E + PV

Ex.3：温度圧力一定	

!

!

E " µN
Ex.4：グランドカノニカル	!
(l )T and "のランダムウォーク	

!

!

WHAM	
  
でSTの	
  
重みを決定	

M 0 M1!

n(E 0 ,V1,!,VL ) = M 0" "! " nm 0 =1m1 =1

exp(# f m 0 ,m1 ,!m L ) =!

(E 0 ,V1,!,VL )m L =1m 0 ,m1 ,!m L

exp( f m 0 ,m1 ,!m L # $ m 0 E % m )m L =1" n(E ,V!,!,V )exp(#$ E )0E 0 ,V1 !VL!

m 0 ,m1 ,!m L

1L

x = {S1,S2 !}x = {q1 ,q2 !,V}x = {q1 ,q2 !,N}E and Vl!0

ML

" "! " Nm 0 =1m1 =1
M1
ML

!

x = {q1 ,q2 !}m&mのランダムウォーク	  	
  計算機のCPUの数（資源）に	
  より使い分けられる。ST、	  STREMはREMより通信が少なくてすむ（レプリカ数すくないから）	

水中のproteinGの折り畳みシミュレーション	
Collaborate	
  with	  Dr.	  T.	
  Yoda	  (Nagahama	  Bio.	
  Univ.),	
  Dr.	
  Yuji	  Sugita	
  (Riken),	
  	
  
	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  Dr.	  T.	
  Nishikawa	
  (Tokyo	
  Tech.	
  Univ.),	
  and	
  Prof.	
  Y.	
  Okamoto	
  (Nagoya	
  Univ.)	
  ••半径50Åの水球中（水分子17,063個）の小タンパク質 protein	
  G（アミノ酸数56,全原子数52,044個）を拡張アンサンブルシミュレーションにより折り畳む	  ProgeinG	
  – Tf=343K(70℃)の2状態フォールディング	  – t	
  (folding)	
  <	
  2	
  msec	  under	  naive	  condiions	
  (Alexander	
  et	
  al.,	
  Biochemistry	
  31,	
  7243	
  (1992))	
  – t	
  (folding)	
  ~	
  5.2	
  msec	  and	  t	
  (iniial	
  collapse)	
  <	
  1msec	
  with	  D/H	
  exchange	
  at	
  pH=4.0,	
  T=278	
  K	
  in	
  the	  presence	  of	
  0.25	
  M	
  GuHCl	
  (Kuszewski	
  et	
  al.,	
  Protein	  Sci.	
  3,	
  1945	
  (1994))	
  – 初期段階（Trp42の埋め込み、C末のベータヘアピンの形成、N末端とへリックスの形成）	
  ••地球シミュレータと東大のT2Kを利用	
  レプリカ交換分子動力学シミュレーション	
  
–
–
–

•

300Kから700Kまでの224コのレプリカ(112node(112*8=896cpu)	
  of	
  Earth	  Simulator)	
  0.1ps毎にレプリカ交換	
  
各6ns	
  

レプリカ交換マルチカノニカルシミュレーション	
  ––WHAMとレプリカ交換シミュレーションの結果から重みを決定	
  56コのレプリカ(224cpu	  of	
  T2K	  in	  東大）	
  ––1ps毎にレプリカ交換	  各13ns	
  初期構造	レプリカ交換シミュレーション	• 1replicaあたり5ns(5*224=1.12μs)*100~1000	
  
	
  (t	
  (folding)	
  〜	
  5.2	
  msec	  	
  and	  	
  t	  (iniial	
  collapse)	
  <	
  1	
  mse)	
  -110000

20
40
60
80
100
120
140
160
180
200
220

-120000-130000

E-140000-150000-160000-170000-180000
0

1e+06

2e+06

3e+06

4e+06
step

5e+06

6e+06

7e+06

8e+06

一番小さいRMSDの構造得られた構造	
0.8
0.6
0.4
0.2
0
1

5 10 15 20 25 30 35 40 45 50 55

55 1
50
45 0.840
35 0.630
25 0.420
15
10 0.2
5
0

天然構造	
0.8
0.6
0.4
0.2
0
1

5 10 15 20 25 30 35 40 45 50 55

55 1
50
45 0.840
35 0.630
25 0.420
15
10 0.2
5
0

ターンの位置やバックボーンの配置が	
  似ている構造が得られた。	
  ヘリックスはできていない。	
  数倍の計算は必要	専用計算機を用いた計算	••••

近年、専用計算機を用いたシミュレーションが盛んに行われている。	
  GPGPU	  (NAMD	  (Shulten),	
  AMBER	  (Pande)	
  ANTON	
  (D.	  E.	
  Shaw)	
  MD-­‐GRAPE3	
  （AMBER(理研)）	
  •

MDGRAPE3を用いて水中のEnkephalinのレプリカ交換MDを行なった。レプリカ交換の部分はperlを用いた。	
  5残基からなるEnkephalinに２４Aの球（1710個の水）	
  Amber94	
  孤立系、Langevin熱浴、36	
  個のレプリカ（275Kから427.1K)	
  各レプリカ4ns	
  (全部で72ns)	
  

•
•
•
•

主成分解析を用いた構造分類	

40

REMD
CTMD

30
20

PCA2

10
0
-10
-20
-30
-40
-40

３つのクラスター	

-30

-20

-10

0
PCA1

10

20

30

緑：2nsの300KのMD	
  赤：REMDの300Kの構造	40時間の比較	• MDGRAPE3(1core)（4ns) -­‐>12時間	
  200ps	  (50A)	
  	
  	
 これまでの計算	
  
Ø
34A

26時間	

40A 74時間		
  (4core) -­‐>12時間	
  80ps	  (Earth)	
  50A 260時間	
	
  	
  	
  (2core)	
  	  -­‐>12時間	
  20ps	  (T2K)	
  
	
  	
  	
  並列化(896core	
  224replicas	
  (Earth)	
  	
  112core	
  56replicas	
  (T2K))	
  	
  1coreあたりは20倍ぐらい高速	
  
GPU？	
  
並列計算機、高速計算機を用いて継続	  	マルチカノニカル分子動力学法による	
  分子クラスターの転移現象の研究	金子敏宏（慶應大）（D３）"秋元琢磨（慶應大）、泰岡顕治（慶應大）、光武亜代理（慶應大）、Xiao Cheng Zeng(U.of Nebraska-Lincoln)"
T. Kaneko, T. Akimoto, K. Yasuoka, A. Mitsutake, and X.C. Zeng, J. Chem. Theo. And Comput. In press	

•••

上空や宇宙空間など低温•低密度条件下では、初期条件に依存して数分子が集合して存在 分子クラスターの性質のサイズ依存性が重要	
  水クラスターの構造が壊れることに対応する比熱のピークが存在	
  水クラスターは一度固相のようになって	
  時間を経ると別の固相に変化する（動的共存現象）	
  水クラスターの転移のサイズ依存性を明らかにする	
  
	
  
	
  
マルチカノニカル分子動力学法を用いて	
  融点を決定	  	
  カノニカル分子動力学法を用いて動的共存の	
  様子を詳しくしらべる	

N=8,12	
  ポテンシャルエネルギーにギャップが存在	  固相のような状態と液相のような状態が時間的に共存	  動的共存において長時間相関が存在	  	
  N=7,11	  両者の区別がつきにくい	
  動的共存において長時間相関が存在しない	  	
  原因ー＞絶対零度でのポテンシャルエネルギーの安定性が影響	スリット型細孔に閉じ込められた	  単純液体の固液相転移	
金子敏宏（慶應大）（D３）"泰岡顕治（慶應大）、光武亜代理（慶應大）、Xiao Cheng Zeng(U.of Nebraska-Lincoln)"	目的とこれまでの研究	
化学物理に基づく全原子の計算機シミュレーションによって蛋白質のフォール
ディングや会合系などの機構について理解する	
フォールディング、会合系のシミュレーション	
  （たくさんの極少値からなる自由エネルギー曲面）	
  • １つの極少状態にトラップされる	
  -­‐>	
  効率の良いサンプリング方法が必要	  – 拡張アンサンブル法をペプチド系へ適用し有効性を調べた	
  – 複合型拡張アンサンブル法の開発を行った	  （近年、多変数版（会合系シミュレーションのため））	
  – 陽に水を取り入れたProteinGの大規模なフォールディングシミュレーション	  – 分子の有限サイズ系への拡張アンサンブル法の適用	  •正確なエネルギー関数が必要（特に、溶媒を正確に取り入れる）	
  – 陽に水分子を加える	  • 近年、専用ボードMD-­‐GRAPE3やGPGPUを用いてシミュレーションを行っている	  – 統計力学に基づく液体論であるRISM	
  理論を用いる	
  • 1D	
  RISM理論を用いた拡張アンサンブルシミュレーション	
  • ３D	
  RISM理論を用いてペプチドと蛋白質の会合系への適用	
  緩和モード解析手法を蛋白質系へ適用	
  – ペプチド系で有効性を調べた	
  – リゾチウムという蛋白質系に適用	
  3次元RISM理論を用いたシミュレーション	
分子研平田グループと共同研究	
•••

統計力学の液体論に基づいた理論	

水などの分子性液体に適用	

溶質の周りの溶媒の3次元密度分布を得ることができる	

3D-RISM Eq.	

h!uv (k ) = % c"uv (k )#"!uv (k ) + % c"uv (k )$ v h"!vv (k )"Closure	"

単原子系では	  Ornstein-­‐Zernike	
  Eqに対応	

h (r) = exp["u!uv (r) / kT + h!uv (r) " c!uv (r)] "1
uv
!

u(r)	
  	
  ：ポテンシャル	  h(r) :	
  the	
  total	
  correlaion	  funcion	  c(r) ：the	  direct	
  correlaion	  funcion	
  
	
  !
	
  	
  	
  	
   ：溶媒の密度	  ! ：The	  intramolecular	
  correlaion	
  matrix	
  2体分布関数はg(r)=h(r) + 1より得られる	
nTINKERプログラムのモンテカルロ部分を改良し,	
  共同研究のもとRISMの部分をエネルギー関数に加える	  n計算時間はかかるが,	
  水の長時間平均分布を得ることができる	  n物理の理論により水の分布を予測できる	• ペプチドと蛋白質の隙間に水の分布	緩和モード解析	
H.	
  Takano	
  and	  S.	
  Miyashita,	
  J.	
  Phys.	
  Soc.	
  J,	
  64	
  (1995)	
  3688.	
  S.	
  Koseki,	
  H.	
  Hirao,	
  and	  H.	
  Takano,	
  J.	
  Phys.	
  Soc.	
  J,	
  66	
  (81997)1631.	
  高分子系で動的性質を調べるのに開発された	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  
-­‐>	
  	
  蛋白質系へ	
  • 構造揺らぎの動的性質を抜き出す。	
  • (PCAは構造揺らぎの静的性質を抜き出す。）	
  	
  ー＞緩和モード解析は主成分解析の動的拡張	  • 緩和モード解析では、緩和率、緩和モードを求めることができる。	
  •X p (t)X q (0) eq = " p,q exp(# $ p t)
X p :p番目に遅い緩和モード	
 " p :	
 p番目に遅い緩和率	
•

時間相関関数を指数的な緩和項で展開する。	
  

!
!
•

!

Ri (t)R j (0) eq = " gi g j exp(#$ p t)!n

遅い緩和モードから得ることができる。	
  主成分解析と緩和モード解析の比較	
• 主成分解析 （平均構造からの構造揺らぎの大きい順）	
  多変数の統計データをある軸に射影したとき、射影したデータの	  分散が最大となる軸を順次決めて行く	  相対座標に関する分散共分散行列の固有値問題に帰着	

C(0)v p = " 2p v pC(0) i, j =< Ri R j > Ri = ri " < ri >" 2p =< (R • v p ) 2 >Ø

緩和モード解析 （構造揺らぎの緩和が遅い順）	
  
N

!
t
X p (Q) = " f pi Ri ( 0 ;Q) Ri (t;Q) :状態Qから出発して時刻t後のRiの期待値	
2
!
i=1

!

時間間隔 t 0 + " および t 0 に対する相対座標に関する平衡時間相関関数の	  一般化固有値問題に帰着	

C(t 0 + " )f p!= exp(#$ p " )C(t 0 )f p!

!

C(t) i, j =< Ri (t)R j (0) >!
緩和モード解析は主成分解析の動的拡張	

シミュレーションと解析手法の詳細	
•系	
  いくつかの準安定構造をもつ真空中のMet-­‐Enkephalin	
  	
  
	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  (Tyr-­‐Gly-­‐Gly-­‐Phe-­‐Met	
  の5残基ペプチド）	
  複数の水素結合のパターンにより構造分類されている （AとB）	
  • MCシミュレーション	
  力場:	
  ECEPP/2、温度:	
  298.15	
  K、２面角を変化、	
  ステップ数:	
  8x１０８MC sweeps、	
  100	  MC	  sweeps毎に構造を保存	
  • 解析について	  各残基の主鎖のCαとCを取り出し、並進と回転の自由度を除く	
  
	
  	
  	
  	
  	
  	
  	
  	
  t0=8000	
  MC	
  sweep	
  tau=2000	
  MC	
  sweep	
  緩和モード解析（RMA）	
  A.Mitsutake,	
  T.	
  Iijima,	
  H.	
  Takano	
  submi?ed	
  to	
  JCP）	

•••高分子の分野で開発された（主成分解析(PCA)を拡張したものに対応）	
  蛋白質系に用いる	  	
  極少状態の分かっているぺプチド系で自由エネルギー曲面を調べる	  – 違う自由エネルギー曲面が得られた	
  – PCAは構造の類似性を調べることができた	
  – RMAは緩和の遅い動きを抜き出すことができた	

PCA	
  

RMA	
  

ー＞より大きな蛋白質系へ適用	
  溶媒の動的動きの解析（実験との比較）	

リゾチームの緩和モード解析 （長井君）	
•
•
•
•
•
•
•
•
•

リゾチームの300Kでの等温シミュレーション	
  AMBER8＋MDGRAPE3を用いる	
  力場	
  ﬀ03	
  原子数：24475原子（リゾチーム:2021原子、Cl-­‐:8原子、水分子:7482	
  分
子）	
  
初期構造：1REX.pdb	
  系：孤立系 water	
  cap=40A	  Langevin熱浴法	  時間刻み：2fs	
  (SHAKE法によりHを固定）	
  全シミュレーション時間：100ns	
  PCAの結果	

RMAの結果	

今後の計画	
• 蛋白質系への緩和モード解析• 拡張アンサンブル法を用いた	
  小蛋白質のフォールディングシミュレーション	  • 会合系のシミュレーション	  – RISM理論を用いて溶媒和を調べる	
  – 会合系の拡張アンサンブル法を既存のソフトに導入	
  
	
  
	
  
ー＞ GPUやMDGRAPE3を用いる	
  	
  • 生体高分子系の計算	

	
  

謝辞	
◆ 拡張アンサンブル法	
  
•
•
•
•
•
•

岡本祐幸教授（名古屋大学）	
  杉田有治準主任研究員（理研、元分子研助手）	
  	
  	
  
Ulrich	  H.E.	
  HANSMANN（ミシガン工大、元分子研助手）	
  依田講師（長浜バイオ大）	
  西川武志（計算科学振興財団）	
  Charles	  L.	
  Brooks	
  III	  (ミシガン大学）	
  ◆ 溶媒効果	
  
•
•
•
•
•

平田文男教授（分子研・総研大）�
丸山博士（分子研）
木下正弘教授（京大・エネルギー理工研）�
今井隆志講師（立命館大）�
原野雄一講師（大阪大）�

◆ 高分子系の解析手法"•

高野宏教授（慶應大）、飯島君（慶應大）、安藤君（慶應大）、長井君（慶應大）"u 粒子系の拡張アンサンブルシミュレーション"•

泰岡顕治（慶應大）、金子敏宏（慶應大）、Xiao Cheng Zeng(U.of Nebraska-Lincoln)"
秋元琢磨（慶應大）"◆ 計算機"
分子科学研究所計算機センター、地球シミュレーター、NAREGI、PETAプロジェクト、T2K(東大）""