特集「生体高分子の揺らぎとダイナミクス—シミュレーションと実験の統計解析—」［研究詳解］

統計数理（2014）
第 62 巻 第 2 号 221–241
c 2014 統計数理研究所


生体高分子系の緩和モード解析
高野 宏†（受付 2013 年 12 月 31 日；改訂 2014 年 7 月 15 日；採択 7 月 16 日）

要旨

詳細釣り合いの成り立つシミュレーション法に対して，時間発展演算子の固有値，固有関数から，扱っている系の緩和率，緩和モードをそれぞれ定義することができる．計算機シミュレー
ションの時系列データから緩和率分布と緩和モードを変分問題に基づいて評価する方法である
緩和モード解析の方法について解説する．さらに，ペプチド系のモンテカルロシミュレーショ
ン，タンパク質系の分子動力学シミュレーションへ応用した例について紹介する．
キーワード：緩和モード，緩和率，緩和モード解析，シミュレーション，時間相関関数．
1.

はじめに

緩和モード解析の方法（Takano and Miyashita, 1995; Koseki et al., 1997; Hirao et al., 1997）
は，計算機シミュレーションで得られた時系列データから，扱っている系の緩和モードと緩和率分布を評価する方法である．従来，ランダムスピン系，高分子系に用いられてきた．本稿では，この方法および生体高分子系への応用例（Mitsutake et al., 2011; Nagai et al., 2013）について説明する．
生体高分子であるタンパク質の機能は，その構造の揺らぎと関連があることが知られている．
タンパク質のシミュレーションで得られたデータから，構造の揺らぎを解析する方法として，従来，主成分分析の方法が広く用いられている．この方法は，構造揺らぎの分散共分散行列の対角化により，分散の大きい揺らぎのモードを順次求めて行く方法で，構造のデータが時系列としてどのように並んでいるかという情報を用いない静的解析である．これに対し，緩和モード
解析の方法は，二つの異なる時間間隔に対する構造揺らぎの時間相関行列に関する一般化固有値問題を解く方法で，時系列としての情報を用いた動的解析である．生体高分子のシミュレー
ションの同じデータに対して，両方法を適用した結果がどのように異なるかは興味深い問題である．
以下では，第 2 節で系の時間発展演算子の固有値，固有関数からどのように緩和率，緩和モードが定義されるかについて説明する．第 3 節では，変分問題に基づいてシミュレーションのデー
タから緩和モード，緩和率を評価する方法である緩和モード解析の方法について説明する．第4 節では，ペプチド系のモンテカルロ・シミュレーションに緩和モード解析を適用した結果について報告し，主成分分析の結果との比較を行う．第 5 節では，タンパク質系の分子動力学シ
ミュレーションに緩和モード解析を適用した例について報告する．この場合，主成分分析で得られた，揺らぎの大きいいくつかの主成分に基づいて緩和モード解析を行う．その際，速い緩和モードの影響を取り除くために行った工夫についても説明する．第 6 節でそれまでの内容を†

慶應義塾大学 理工学部：〒 223–8522 神奈川県横浜市港北区日吉 3–14–1

統計数理

222

第 62 巻

第 2 号 2014

まとめる．緩和モード解析の方法では重要な物理量の線形結合を変分問題の試行関数としている．固有値方程式と変分問題の関係について付録 A で，試行関数の選び方について付録 B で説明する．
2.

緩和モードと緩和率

詳細釣り合いの条件を満たす三つのシミュレーション方法に対して，時間発展演算子の固有値，固有関数から緩和率，緩和モードがどのように定義されるかについて説明する．
最初のシミュレーション方法は，詳細釣り合いの条件を満たした連続時間のモンテカルロ法である．この方法では，時刻 t に系が状態 Q をとる確率 P (Q; t) の時間発展がマスター方程式

∂
(2.1)
Γ (Q|Q )P (Q ; t)
P (Q; t) = −∂tQで与えられる．ここで，Γ (Q|Q ) は時間発展演算子（遷移行列）であり，平衡分布関数 Peq (Q) に対し詳細釣り合いの条件(2.2)

Γ (Q|Q )Peq (Q ) = Γ (Q |Q)Peq (Q)を満たすように選ばれる．
次のシミュレーション方法は，Brown 動力学法で，N 個の粒子からなる系において，粒子の位置を ri (i = 1, . . . , N ) で表し，粒子間の相互作用のポテンシャルを U ({rj }) とするとき，粒子の位置の時間発展は Langevin 方程式


d ri
1
∂U ({rj })
+ wi
=
−
(2.3)
dt
ζ
∂ri
で与えられる．ここで w i は i 番目の粒子に働くランダム力で，時刻 t におけるその α 成分を
wiα と表すと，それらは(2.4)

wiα (t) = 0

および
(2.5)

wiα (t)wjβ (t ) = 2ζkB T δi,j , δα,β δ(t − t )を満たす白色 Gauss 過程である．ここで，ζ は抵抗係数で，速度に比例した抵抗力の比例定数である．kB は Boltzmann 定数，T は系の絶対温度を表す．このとき，時刻 t に系が状態
Q = t (t r 1 , . . . , t r N ) をとる確率 P (Q; t) の時間発展は Smoluchowski 方程式
(2.6)

∂P
= −Γ P
∂t


 ∂
1
∂
∂U
·
+
=
kB T
P
∂r i ζ
∂r i∂r iiで与えられる．時間発展演算子 Γ は平衡分布 Peq (Q) ∝ exp[−U/(kB T )] に対し詳細釣り合いの
条件
(2.7)

Peq (Q )Γ (Q)δ(Q − Q ) = Peq (Q)Γ † (Q)δ(Q − Q )を満たす．ここで Γ (Q)δ(Q − Q ) と Γ † (Q)δ(Q − Q ) は Γ とその共役 Γ † が δ(Q − Q ) の Q にのみ作用することを表す．Γ (Q)δ(Q − Q ) = Γ (Q|Q )，Γ † (Q)δ(Q − Q ) = Γ (Q |Q) と表す演算子の行列表現を用いれば，この詳細釣り合いの条件は式（2.2）と同じである．

生体高分子系の緩和モード解析

223

三つ目のシミュレーション方法は，Langevin 熱浴を温度制御に用いた分子動力学法で，簡単のため，同じ質量 m，同じ抵抗係数 ζ をもった N 粒子系を用いて説明する．i 番目の粒子の位
置 r i と速度 v i の従う Langevin 方程式は
m

(2.8)

d vi
∂U ({r j })+ wi,
= −ζv i −
dt
∂r i

d ri
= vi
dt

で与えられる．ここで，U ({r j }) は粒子間の相互作用のポテンシャルである．w i は i 番目の粒子に働くランダム力で，Brown 動力学法で用いたものと同じ白色 Gauss 過程である．このとき，
時刻 t に系が状態 Q = t (t r 1 , . . . , t r N , t v 1 , . . . , t v N ) をとる確率 P (Q; t) の時間発展は Kramers
方程式
(2.9)

∂P
= −Γ P
∂t


 ∂
∂
1 ∂V
∂
ζ
kB T ∂· vi +
·
+
·
=
vi +
−
P
∂r i
∂v i m ∂r i
∂v i mm ∂v iiで与えられる．時間発展演算子 Γ は平衡分布 Peq (Q) ∝ exp[− kB1T ( m2
詳細釣り合いの条件

2i v i + U ({r j }))] に対しPeq (Q )Γ (Q)δ(Q − Q ) = Peq (Q)Γ † (Q)δ(Q − Q )(2.10)

を満たす（Risken, 1996）．ここで Q = t (t r 1 , . . . , t r N , −t v 1 , . . . , −t v N ) は状態 Q の速度を反転した状態である．行列表現を用いれば，この詳細釣り合いの条件は
Γ (Q|Q )Peq (Q ) = Γ (Q |Q)Peq (Q)

(2.11)と書ける．ここで，平衡分布は
Peq (Q) = Peq (Q)(2.12)

を満たす．
以上の時間発展演算子の固有関数としての緩和モードについて説明する．時刻 t に系が状態
Q をとる確率 P (Q; t) の時間発展の方程式を行列表現で表すと（2.1）式となる．モンテカルロ法，Brown 動力学法では，速度の自由度がないため，Q = Q となるので，分子動力学法の場合とで与えられる．時間発展演算子 Γ (Q|Q ) の固有値問題
合わせて，詳細釣り合いの条件は（2.11）を考える．固有値 λn に対する左固有関数を φn (Q)，右固有関数を ψn (Q) とする．即ち，
(2.13)φn (Q)Γ (Q|Q ) = λn φn (Q ),Q

(2.14)

Γ (Q|Q )ψn (Q ) = λn ψn (Q).

Q

右固有関数を ψn (Q) = φ̂n (Q)Peq (Q) と書くとき，φ̂n (Q) は左固有関数 φn (Q) を使って φ̂n (Q) =（2.13）式の両φn (Q) と書ける．このことは次のようにしてわかる．左固有関数の固有値方程式

辺に Peq (Q ) をかけて
(2.15)φn (Q)Γ (Q|Q )Peq (Q ) = λn φn (Q )Peq (Q )Qを得る．この左辺は詳細釣り合いの条件（2.11）式を使って，

(2.16)φn (Q)Γ (Q|Q )Peq (Q ) =φn (Q)Γ (Q |Q)Peq (Q)
Q

Q

統計数理

224

=



第 62 巻Γ (Q |Q)φn ((Q))Peq (Q) =第 2 号 2014Γ (Q |Q )φn (Q )Peq (Q )

Q

Q

と変形できる．一方（2.15）
式の右辺は (Q) = Q と（2.12）式を用いてλn φn (Q )Peq (Q ) = λn φn ((Q ))Peq (Q )(2.17)と変形できる．この式の右辺と（2.16）
式の最後の式を等号で結び，Q を Q と，Q を Q と書き換えて(2.18)
Γ (Q|Q )φn (Q )Peq (Q ) = λn φn (Q)Peq (Q)Q

を得る．これは右固有関数の固有値方程式（2.14）式で，ψn (Q) = φn (Q)Peq (Q) であることを示す．即ち，
ψn (Q) = φ̂n (Q)Peq (Q) = φn (Q)Peq (Q).

(2.19)

固有関数は規格直交関係

(2.20)φm (Q)ψn (Q) =φm (Q)φ̂n (Q)Peq (Q) = φm φ̂n  = δm,nQ

Q

を満たすようにとる．ここで，A は物理量 A の平衡状態での平均を表す．このとき，時間間隔 t だけ離れた φm (Q) と φ̂n (Q) の時間相関関数は
(2.21)
φm (t)φ̂n (0) =φm (Q)Tt (Q|Q )φ̂n (Q )Peq (Q )
=


Q

=



Q

Q

φm (Q)e−Γ t (Q|Q )φ̂n (Q )Peq (Q ) =

Q


Q

φm (Q)e−Γ t (Q|Q )ψn (Q )Q

φm (Q)e−λn t ψn (Q) = δm,n e−λn tQと与えられる．ここで，Tτ (Q|Q ) = e−Γ τ (Q|Q )(2.22)

は時刻 t で系が状態 Q にあるとき時刻 t + τ で系が状態 Q にある条件付き確率である．
2 つの物理量 A(Q) と B(Q) が固有関数 φn (Q) と φ̂n (Q) を使って
(2.23)
A(Q) =
an φn (Q), B(Q) =b̂n φ̂n (Q)nnのように固有関数展開できるとき，A(Q) と B(Q) の時間相関関数は
(2.24)
A(t)B(0) =an b̂n exp(−λn t)nで与えられる．この性質から，固有関数 φn (Q) と φ̂n (Q) を緩和モード，固有値 λn をその緩和率と呼ぶことができる．
3.

緩和モード解析の方法

緩和モード解析の方法は，時間発展演算子に対する固有値問題と等価な変分問題を考え，適切な試行関数を選ぶことで，緩和モード，緩和率を評価する方法である．まず，時間発展演算

生体高分子系の緩和モード解析

225

子そのものではなく，条件付き確率（2.22）に対する固有値方程式
(3.1)
φn (Q)Tτ (Q|Q ) = e−λn τ φn (Q ),
Q



(3.2)

Tτ (Q|Q )ψn (Q ) = e−λn τ ψn (Q)Q

を考える．この固有値方程式と等価な変分問題は
Q

R=

(3.3)



Q φn (Q)Tτ (Q|Q )ψn (Q )Q φn (Q)ψn (Q)=φn (t)φ̂n (0)
φn φ̂n 

に対し，δR = 0 で与えられる（Talkner, 1994; Hirao et al., 1997; Drozdov, 2000）．付録 A で固有値方程式と変分問題について説明した．δR = 0 となる φn , φ̂n が固有関数で，そのときの Rの値が exp(−λn τ ) となる．この変分問題に対し，重要と思われる物理量の線形結合を試行関数として用いることで，計算機シミュレーションのデータから緩和モードと緩和率を評価することができる．
緩和モードの計算法の例として，N 個のセグメントからなる孤立線状高分子のシミュレーショ
ンを考える．i 番目のセグメントの位置を r i ，速度を v i とし，その重心からの相対位置を Ri ，相対速度を RN +i とする．即ち，(3.4)

R = t (t R1 , . . . , t RN , t RN +1 , . . . , t R2N )= t (t r 1 − t r c , . . . , t r N − t r c , t v 1 − t v c , . . . , t v N − t v c ).

ここで r c = N1 i r i と v c = N1 i v i は，それぞれ重心の位置と速度である．φn (Q)，φ̂n (Q)に対する試行関数（近似的固有関数）
として，それぞれ，
(3.5)
X p (Q) =fp,i Ri (t0 /2; Q), X̂ p (Q) =fp,i Ri (t0 /2; Q)iiを考える．ここで，Ri (t; Q) =

(3.6)



Ri (Q )Tt (Q |Q)Q

は時刻 0 で系が状態 Q にあるとき時刻 t での Ri の期待値である．パラメター t0 は，Ri の中に含まれる緩和の速いモードの寄与を相対的に小さくするためのものである．付録 B で t0 の役
割について説明した．fp,i は変分パラメターである．この試行関数は，左右固有関数の間の関係（2.19）式 X̂ p (Q) = X p (Q) を満たしている．
式に代入すると，fp,i に対する変分問題 δR = 0 は一般化固有値問題
これらの試行関数を
（3.3）


(3.7)
fp,i Ci,j (t0 + τ ) = e−λp τfp,i Ci,j (t0 ),
(3.8)



i

i

Ci,j (t0 + τ )j fp,j = e

j

−λp τ



Ci,j (t0 )j fp,jjと等価になる．ここで，
(3.9)

i =

1,
−1 ,

i = 1, . . . , N
i = N + 1, . . . , 2Nである．Ci,j (t) = 13 Ri (t) · Rj (0) は時間相関行列であり，詳細釣り合いの条件（2.11）式より，（3.9）式で定義された i を用いて Ci,j (t) = i j Cj,i (t) という対称性を持っている．この対称性

統計数理

226

第 62 巻

第 2 号 2014より，fp,i に関する一般化固有値問題（3.7）
式と i fp,i に関する一般化固有値問題（3.8）式が等価であることがわかる．
式は X p (Q)，X̂ p (Q) の α 成分を，それぞれ Xp,α (Q)，X̂p,α (Q)
固有関数間の規格直交関係（2.20）と書くと
Xp,α X̂q,β  = δp,q δα,β(3.10)

で与えられるが，これは fp,i に対する条件
(3.11)
fp,i Ci,j (t0 ) j fq,j = δp,qijに相当する．このとき，近似的緩和モードの時間相関関数として
Xp,α (t)X̂q,β (0)  δp,q δα,β exp(−λp t)(3.12)

が成立すると考えられる．
式より，Ri (t0 /2; Q) の線形結合で試行関数 X p (Q) と X̂ p (Q) を与える
規格直交関係（3.11）（3.5）式の逆変換は

(3.13)Ri (t0 /2; Q) =gi,p X p (Q), i Ri (t0 /2; Q) =ĝi,p X̂ p (Q)ppとなる．ここで，gi,p , ĝi,p は，それぞれ
(3.14)
gi,p =
Ci,j (t0 ) j fp,j ,ĝi,p =j

fp,j Cj,i (t0 ) = i gi,pjで与えられる．これより Ri の相関行列は，t ≥ t0 のとき，1
1 
Ci,j (t) = Ri (t) · Rj (0) =Ri (Q) · Tt (Q|Q )Rj (Q )Peq (Q )
3
3 Q 
Q

1 Ri (Q)Tt0 /2 (Q|Q ) · Tt−t0 (Q |Q )Tt0 /2 (Q |Q )Rj (Q )Peq (Q )=3 Q   
Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

1 
Ri (t0 /2; Q ) · Tt−t0 (Q |Q )Tt0 /2 (Q |Q )Rj (Q )Peq (Q )=3   

1 Ri (t0 /2; Q ) · Tt−t0 (Q |Q ){j Rj (Q )Tt0 /2 (Q |Q )}Peq (Q )=3   
1 Ri (t0 /2; Q ) · Tt−t0 (Q |Q )j Rj (t0 /2; Q )Peq (Q )=3   
1 
gi,p X p (Q ) · Tt−t0 (Q |Q )ĝj,q X̂ q (Q )Peq (Q )=3    pq

1
gi,p ĝj,q X p (t − t0 ) · X̂ q (0)
3
p
q



gi,p ĝj,p exp[−λp (t − t0 )] =g̃i,p g̃ˆj,p exp(−λp t)

=

p

p

のように与えられる．ここで，(3.15)
g̃i,p = eλp t0 /2
Ci,j (t0 ) fˆp,j ,jg̃ˆi,p = eλp t0 /2jfp,j Cj,i (t0 ) = i g̃i,p生体高分子系の緩和モード解析

227である．途中の変形で詳細釣り合いの条件（2.11）式から得られる結果
T (Q|Q )Peq (Q ) = T (Q |Q)Peq (Q)

(3.16)

および（2.12）式を用いている．（3.15）
式の最後の行は，Ri (Q) の（2.23）式のような展開がg̃ˆi,p X̂ p (Q)(3.17)Ri (Q) g̃i,p X p (Q) ppのように与えられると見なすことができる．
以上，相関行列 Ci,j (t0 ) と Ci,j (t0 + τ ) をシミュレーションで計算し一般化固有値問題（3.7）
式と緩和率 λp を評価することができる．その結果
式，（3.8）式を解くことで，緩和モード（3.5）に基づいて，（3.15）
式を用いて t ≥ t0 に対する時間相関行列を計算し，実際に計算機シミュレー
ションで求めた結果とくらべることで，得られた緩和モード，緩和率が妥当なものであるかどうかを検証することができる．
緩和モード解析の結果の例として，Langevin 熱浴を温度制御に用いた孤立線状高分子の分子
動力学シミュレーションへ適用した Hirao et al.（1997）の結果を紹介する．
孤立線状高分子の N 個のセグメント間の相互作用ポテンシャルは Kremer and Grest（1990）と同じもので，全セグメント間に働く斥力相互作用（Lennard-Jones ポテンシャルの斥力部分）と高分子鎖上で隣接するセグメント間に働く引力相互作用から成る．温度や抵抗係数等のパラメと同じである．このポテンシャルとパラメターでは，セグメン
ターも Kremer and Grest（1990）
ト間の排除体積相互作用が取り入れられることが知られている．流体力学的相互作用を無視して，（2.8）式に基づいて，セグメント数 N = 10, 20, 40 の場合にシミュレーションを行った．平衡
4Lennard-Jones 単位）のデータを用いて時間相関
化後の時間 4 × 10（単位は斥力ポテンシャルの（3.8）式を解き，緩行列 Ci,j (t) を計算し，t0 = 1.6, τ = 0.4 として，一般化固有値問題（3.7）式，
和モード，緩和率を評価した．得られた結果から（3.15）式を用いて重心からの相対位置の自己相（i = 1, . . . , N ）と重心からの相対速度の自己相関関数 Ci,i (t)（i = N + 1, . . . , 2N ）と関関数 Ci,i (t)を計算し，実際にシミュレーションから求めた結果と比較したのが図 1 である．減衰振動的な挙動も良く再現されており，得られた緩和モード，緩和率が妥当であることを示している．分子動力学法の場合，振動的な振舞いも出てくるので，緩和率は複素数となる．その分布を図 2 に示
す．緩和率の実部が小さいときは，緩和率が実数になり，純粋緩和となっていることがわかる．
緩和率が実数となっている緩和率について，緩和率が小さい順に緩和率 λp の番号 p を 1 からつける．このとき，λp ∝ (p/N )x というベキ乗則的振舞いが見られた．フィッティングで得られた指数は x  2.2 となっており，理論の予想値 x  2.18
（Doi and Edwards, 1986）に近い値が得

式の展開係数 g̃i,p の i 依存性は，理論の予想 g̃i,p ∝ cos i − 12 pπ/N
られている．また，（3.17）に良く合っていた．
4.

ペプチド系への応用

生体高分子系への応用例の 1 つとして，真空中の Met-Enkephalin のモンテカルロ・シミュレー
ションへ緩和モード解析を適用した Mitsutake et al.（2011）の結果を紹介する．Met-Enkephalin
は Tyr-Gly-Gly-Phe-Met の 5 残基からなるペプチドで，真空中では複数の準安定状態をもつことが知られている（Mitsutake et al., 1998）．そのため，安定状態，準安定状態間の遷移を見るのに適した系である．シミュレーションの詳細は原論文に譲り，ここでは，主成分分析と緩和モード解析という 2 つの方法の結果を示す．
解析に用いる変数は

228

統計数理

第 62 巻

第 2 号 2014

図 1．N = 10 の場合の Ci,i (t) の t 依存性．実線は緩和モード解析の結果に基づいて（3.15）式を用いて計算した結果を，黒丸はシミュレーションから直接評価した結果を表す．
（a）には重（上から i = 1, 2, 3, 4, 5）
心からの相対位置の自己相関関数 Ci,i (t) = CN −i+1,N −i+1 (t)が，（b）には重心からの相対速度の自己相関関数 Ci,i (t) = C3N −i+1,3N −i+1 (t)（t  8で下から i = 11, 12, 13, 14, 15）が示されている．Reprinted with permission from H.
Hirao, S. Koseki, and H. Takano, J. Phys. Soc. Jpn. 66, 3399（1997）. Copyright
1997, The Physical Society of Japan.

図 2．λp の分布．菱形，正方形，丸がそれぞれ N = 10, 20, 40 の結果を表す．Reprinted withpermission from H. Hirao, S. Koseki, and H. Takano, J. Phys. Soc. Jpn. 66, 3399
（1997）. Copyright 1997, The Physical Society of Japan.

(4.1)

R = t (t R1 , . . . , t RN ) = t (t r 1 − t r 1 , . . . , t r N − t r N )

である．ここで，r i は i 番目の原子の重心からの相対位置座標であり，r i  はその平均位置である．平均位置で作られる構造を平均構造と呼ぶことにする．シミュレーションから得られた多数の構造のサンプルのそれぞれに対し，1 つの構造と平均構造の 2 乗変位が最小になるように，その構造を並進，回転するという操作をしたものを用いる．さらに，操作後の多数の構造を平均したものが平均構造と一致するように平均構造を選ぶ．扱う原子は，主鎖上の 10 個の Cα , C

生体高分子系の緩和モード解析

229

原子とする．即ち N = 10 である．
について説明する．主成分分析では分散共分散行列
まず，主成分分析（PCA）
C = R tR

(4.2)

に対する固有値問題
CF p = Λp F p(4.3)

を規格直交条件
t(4.4)

F p F q = δp,qのもとで解く．p 番目に大きい固有値を Λp とする．固有ベクトル F p は第 p 主成分軸と呼ばれ，この方向に R(Q) を射影した
Φp (Q) = t F p R(Q)(4.5)

は第 p 主成分と呼ばれる．Λp は第 p 主成分の分散を与える．また，R(Q) は F p を使って

(4.6)
R(Q) =
Φp (Q) F p
pのように展開できる．
は時間相関行列
これに対し，緩和モード解析（RMA）
C(t) = R(t) tR(0)

(4.7)

に対する一般化固有値問題(4.8)

C(t0 + τ )f p = e−λp τ C(t0 )f p

を規格直交条件t

f p C(t0 )f q = δp,q(4.9)

のもとで解く．今はモンテカルロ・シミュレーションの場合を扱っているので，φn (Q) = φ̂n (Q)である．p 番目に小さい緩和率を λp とする．対応する緩和モード Xp (Q) と R(Q) の第 i 成分
Ri (Q) の関係は

(4.10)Ri (Q) g̃i,p Xp (Q)pで与えられる．これより，R(Q) は g̃i,p を第 i 成分とするベクトル g̃ p を用いて

(4.11)
R(Q) Xp (Q) g̃ p =Yp (Q) g̃ p /|g̃ p |ppと展開できる．ここで，(4.12)Yp (Q) = Xp (Q)|g̃ p |は PCA の Φp (Q) に対応する量である．
PCA で得られた分散 Λp と RMA で得られた緩和率 λp を図 3 に示す．PCA では p = 1 だけが特に大きい分散を与えるのに対し，RMA では p = 1 と 2 が特に小さい緩和率を与えている．
PCA で得られた F p と RMA で得られた g̃ p /|g̃ p | を p = 1, 2, 3 の場合に図 4 に示す．F p の p = 1

230

統計数理

第 62 巻

第 2 号 2014

図 3．PCA で得られた分散 Λp の分布（a）と RMA で得られた緩和率 λp の分布（b）．Λp の
単位はÅ2 ，λp の単位は自由度あたりモンテカルロ・ステップの逆数．Reprinted withpermission from A. Mitsutake, H. Iijima, and H. Takano, J. Chem. Phys. 135,
164102（2011）. Copyright 2011, American Institute of Physics.

図 4．PCA で得た F p と RMA で得た g̃p /|g̃p |．横軸は主鎖上の原子の番号，縦軸はその原子に対応した F p と g̃p /|g̃p | の x, y, z 成分の値．濃い灰色，中位の灰色，薄い灰色が（a）に，p = 2 を（b）に，p = 3 を（c）に，それぞれ x, y, z 成分を表す．F p の p = 1 を
g̃p /|g̃p | の p = 1 を（d）に，p = 2 を（e）に，p = 3 を（f）に示した．

と g̃ p /|g̃ p | の p = 1 が似ていることがわかる．F p の p = 2 と g̃ p /|g̃ p | の p = 3 が似ているが，g̃ p /|g̃ p | の p = 2 に似ているものはない．

生体高分子系の緩和モード解析

231

図 5．PCA で得られた自由エネルギー地形．各グリッドに入る確率を P とするとき，− ln Pの値を図の右に示した濃度で表している．Φ2 対 Φ1 を（a）に，Φ3 対 Φ2 を（b）に示した．図中の A-H は極小点を示す．Reprinted with permission from A. Mitsutake,
H. Iijima, and H. Takano, J. Chem. Phys. 135, 164102（2011）. Copyright 2011,
American Institute of Physics.

図 6．RMA で得られた自由エネルギー地形．各グリッドに入る確率を P とするとき，− ln Pの値を図の右に示した濃度で表している．Y2 対 Y1 を（a）に，Y3 対 Y2 を（b）に示した．図中の A-H は，図 5 中の極小点 A-H がこの図でどこに相当するかを示している．Reprinted with permission from A. Mitsutake, H. Iijima, and H. Takano, J.
Chem. Phys. 135, 164102（2011）. Copyright 2011, American Institute of Physics.

次に自由エネルギー地形について考える．Φp (Q) または Yp (Q) を X ，Φq (Q) または Yq (Q) をY として，シミュレーション中に X と Y のとる値から，その確率密度分布 P (X, Y ) を求め，無次元化した自由エネルギー F (X, Y ) = − ln P (X, Y ) を求め，X-Y 平面上の等高線図を描いたものが自由エネルギー地形である．図 5 に p = 1, 2, 3 の Φp (Q) で調べた自由エネルギー地形を，
図 6 に p = 1, 2, 3 の Yp (Q) で調べた自由エネルギー地形を示す． PCA で得られた自由エネル
ギー地形の極小点を図 5 の中に A-H で示した．これらの極小点が図 6 中のどこに相当するかを（a）
中で，深い極小は A, C, D で，B は A と C の間に位置し，E は
図 6 中に A-H で示した．図 6
C に属し，F，G，H は D の一部となっている．第 1 緩和モード g̃ 1 /|g̃ 1 | は第 1 主成分軸 F 1 と同様，A と C または A と D の間の遷移を表している．第 2 緩和モード g̃ 2 /|g̃ 2 | は C と D の間の

統計数理

232

第 62 巻

第 2 号 2014

図 7．PCA で得られた自由エネルギー極小の位置での Tyr1 の側鎖の配座の比較．B（濃灰色）（淡灰色）の比較を（a）に，C（濃灰色）と D（淡灰色）の比較を（b）に，E（濃灰色）と FとA（淡灰色）の比較を（c）に示した．それぞれの組について，Tyr1 以外の残基の主鎖部分の平均 2 乗変位が最小となるように配置してある．Reprinted with permission from
A. Mitsutake, H. Iijima, and H. Takano, J. Chem. Phys. 135, 164102（2011）.
Copyright 2011, American Institute of Physics.

遷移を表している．同時にこのモードは残基 Tyr1 の側鎖の配座の変化を記述している．図 7 はA と B 中の配座の比較，D と C 中の配座の比較，F と E 中の配座の比較を示している．それぞれの組で，Tyr1 以外の残基の主鎖部分間の平均 2 乗変位が最小となるように配置している．B,
C, E では Tyr1 の側鎖は主鎖の左側に，A, D, F では右側にある．B, C, E は Y2 が正の領域に，A, D, F は Y2 が負の領域に存在している．図 4（e）に描かれた g̃ 2 /|g̃ 2 | で，Tyr1 に対応する第 1成分の振幅が大きいことがわかる．PCA で得られた F p には g̃ 2 /|g̃ 2 | に似たものはなく，PCA
では抽出できなかった遅い運動，即ち，Tyr1 の側鎖の遅い運動を RMA により抽出できたことがわかる．
また，図 5 と図 6 を比較すると，PCA の自由エネルギー地形では，極小点が原点を取り囲むように並ぶのに対し，RMA の自由エネルギー地形では，極小点が軸に沿って並ぶ傾向がある．
これは，安定状態，準安定状態間の遷移を正しく記述できているためと考えられる．
5.

タンパク質系への応用

生体高分子系へのもう 1 つの応用例として，水中のヒトリゾチームの分子動力学シミュレー
ションへ緩和モード解析を適用した Nagai et al.（2013）の結果を紹介する．ヒトリゾチームは
130 残基，2021 原子からなるタンパク質で，細菌の細胞壁にある多糖類の加水分解を行うという機能をもった酵素である．天然構造で 4 個の α ヘリックスと 1 個の逆平行 β シートをもっている．活性部位にはポケットがあり，このポケットの開閉運動で基質を捕捉することが知られている．ヒトリゾチームに 7482 個の水分子と 8 個の Cl− を加えた 24475 個の全原子に対し，Langevin 熱浴を使った等温分子動力学シミュレーションを温度 300 K で 100 ns 実行した．
式の R であるが，この節では，リゾチームの水素以外の解析に用いる変数は前節と同じ（4.1）
全原子 1029 個を対象とする．前節と同様に，全原子の構造の多数のサンプルそれぞれを，平均構造との 2 乗変位が最小になるように並進，回転したものを用いる．
前節同様，主成分分析を行う．分散共分散行列（4.2）式に対して固有値問題（4.3）式を規格直交条件（4.4）のもとで解く．固有値 Λp は大きい順に番号をつける．第 p 主成分 Φp (Q) は（4.5）式で与えられる．
緩和モード解析を前節と同様に行うことを考える．時間相関行列（4.7）式に対して一般化固有式のもとで解く．分子動力学シミュレーションではあるが，値問題（4.8）式を規格直交条件（4.9）

生体高分子系の緩和モード解析

233

位置のみを扱っているので，時間相関行列 C(t) は対称行列となり，前節のモンテカルロ法の場合と同様の解析となる．これは，速度の緩和が位置に比べて非常に速いため，位置の遅い緩和だけに興味があるときは妥当な解析法である．
しかし，上記の方法はそのままでは使えない．それは，緩和モード解析の方法では，時間相関行列が比較的高い統計精度で評価されている必要があるからである．このことは，次のように理解できる．1 成分の場合の緩和モード解析を考えてみる．1 成分の確率過程を r とし，時間相関関数 C(t) = r(t)r(0) に対し，一般化固有値問題 C(t0 + τ )f = exp(−λτ )C(t0 )f を規格化条件
f C(t0 )f = 1 の下で解く．その結果は λ = {− ln C(t0 + τ ) + ln C(t0 )}/τ , f = C(t0 )−1/2 で与えられる．緩和モードは X = f r(t0 /2; r) = C(t0 )−1/2 r(t0 /2; r) で与えられ，g = C(t0 )f = C(t0 )1/2
を用いて逆変換 r(t0 /2; r) = gX = C(t0 )1/2 X が与えられる．これより，r(t0 )r(0) = C(t0 ) および r(t0 + τ )r(0) = C(t0 + τ ) が再現される．また，g̃ = exp(λt0 /2)g = exp(λt0 /2)C(t0 )1/2 により r(t)r(0)  g̃ 2 exp(−λt) が与えられる．この式は t = t0 と t = t0 + τ では C(t0 ) と C(t0 + τ )を正しく与える．λ に関する結果は，C(t) 対 t の片対数プロットで，t = t0 と t = t0 + τ に対する
C(t) の 2 点を通る直線の傾き −λ を求めていることになる．g に関する結果は，傾き −λ の直線
が C(t0 ) を通るように g を決めていることになる．まとめると，この方法は C(t) = g̃ 2 exp(−λt)
を仮定して，ln C(t) 対 t の片対数プロットで，(t0 , ln C(t0 )) と (t0 + τ, ln C(t0 + τ )) の 2 点を通る
ように直線を引き，傾き −λ と振幅 g̃ 2 を決めていることになる．これより，C(t0 ) と C(t0 + τ )の値が比較的高い統計精度で求まっていないと，λ と g̃ 2 が良い精度で求まらないことになる．
多成分の場合も同様で，緩和モード解析の方法は，Ci,j (t) = p g̃i,p g̃j,p exp(−λp t) の形を仮定
して，Ci,j (t0 )，Ci,j (t0 + τ ) の値から λp と g̃i,p を決めることと等価である．従って，Ci,j (t0 )，Ci,j (t0 + τ ) の値が比較的高い統計精度で求まっている必要がある．この要請のため，緩和モード解析ではタンパク質の全原子の巨大な自由度を直接扱うことはできず，何らかの方法で自由度を逓減する必要がある．
自由度逓減の方法として主成分分析の結果を用いる方法である主成分緩和モード解析（PCRMA）の方法について説明する．この方法では，まず全原子に対して主成分分析を行う．揺らぎの大
きい上位の Nc 個の主成分によって張られる自由度に対して緩和モード解析を行う．このとき，試行関数は上位主成分(5.1)

Φc = t (Φ1 , Φ2 , . . . , ΦNc ),と変分パラメター(5.2)

f p = t (fp,1 , fp,2 , . . . , fp,Nc ),

を用いて
(5.3)

Xp (Q) = t f p Φc (t0 /2; Q).

と与えられる．この試行関数より，相関行列(5.4)

C(t) = Φc (t)t Φc (0)に対する一般化固有値問題（4.8）
式を規格直交条件（4.9）式のもとで解くことになる．ここでは，位置のみを扱っているので，f̂ p = f p が成立する．
上記の方法を試みたところ，主成分の自己相関関数をうまく再現できなかった．用いた発展
時間 t0 の値においては，与えた変分パラメターでは表現できない速い緩和モードが，試行関数の中に含まれていたためと考えられる．言い換えると，上位の主成分には多くの速い緩和モードが含まれており，緩和モードの良い近似になっていない．そこで，2 つの発展時間 t1 と t2 を

統計数理

234

第 62 巻

第 2 号 2014

使った緩和モード解析を導入する．この場合，試行関数は変分パラメター
1
1
1
f 1p = t (fp,1, fp,2, . . . , fp,N
),
c

(5.5)

2
2
2
f 2p = t (fp,1, fp,2, . . . , fp,N
),
c

を用いて
†

Xp (Q) = t f 1p Φc (t1 /2; Q) + t f 2p e−Γ t2 /2 Φc (t2 /2; Q)(5.6)と与えられる．このように試行関数をとるのは，Φc (t1 /2; Q) と Φc (t2 /2; Q) では，消去したい
速い緩和モードの振幅が異なっているため，これらの線形結合で速い緩和モードを表現できると考えるからである．先程の 1 成分の例でいえば，単指数関数が良い近似になってない時間領域で，r(t)r(0) = g̃12 exp(−λ1 t) + g̃22 exp(−λ2 t) のように 2 つの指数関数の和の形を仮定し，4つの t の値における C(t) の値から 4 つのパラメター g̃12 , g̃22 , λ1 , λ2 を決めることに相当する．試行関数（5.6）式から一般化固有値問題
C̃(τ )f p = e−λp τ C̃(0)f p(5.7)と，規格直交条件t

f p C̃(0)f q = δp,q(5.8)

が得られる．ここで，t

f p = (tf 1p , tf 2p )(5.9)

であり，C̃(t) は（5.4）
式で定義された C(t) を用いて⎛⎞
2
C(t1 + t)
C( t1 +t
+ t)
2⎠
C̃(t) = ⎝(5.10)
1
+
t)
C(t
+
t)
C( t2 +t22と定義される．逆変換は
g p = t (tg 1p , tg 2p ) = C̃(0)f p(5.11)を用いて(5.12)Φc (t1 /2; Q)Φc (t2 /2; Q)

=

2N
c

g p Xp (Q)p=1

で与えられる．このとき，Φc (Q) の展開としては(5.13)Φc (Q) 2Ncexp(λp t1 )g 1p Xp (Q),p=1

Φc (Q) 2Ncexp(λp t2 )g 2p Xp (Q)

p=1という 2 種類が考えられる．ここでは 2 つをを平均した展開係数(5.14)

1
2g̃ av
p = {exp(λp t1 )g p + exp(λp t2 )g p }/2を用いた展開(5.15)Φc (Q) 

2N
c

g̃ av
p Xp (Q)p=1を使うことにする．この展開により，（5.4）
式で定義された主成分の相関行列は

生体高分子系の緩和モード解析

235

図 8．PCRMA で得られた緩和時間 λ−1p の分布．Reprinted with permission from T. Nagai,
A. Mitsutake, and H. Takano, J. Phys. Soc. Jpn. 82, 023803（2013）. Copyright
2013, The Physical Society of Japan.

(5.16)

C(t) 

2N
c

t av
g̃ av
g̃ p exp(−λp t)pp=1のように再構成される．
（4.6）式のように表される．上位の主成分だけ
位置の揺らぎを表す変数 R は主成分を使ってを使って，小さい揺らぎを無視したという意味で粗視化された位置の揺らぎを表す変数 Rc を
(5.17)Rc (Q) =NcΦp (Q) F p = F c Φc (Q)p=1で定義する．ここで(5.18)

F c = (F 1 , . . . , F Nc )である．（5.17）式に（5.15）
式を代入することにより，Rc は(5.19)

h̃p = F c g̃ av
pを用いて(5.20)

Rc 

2N
c

Xp (Q)h̃pp=1と表される．（4.11）と同様，h̃p /|h̃p | が第 p 緩和モードの運動の方向を表し，Xp (Q)|h̃p | が PCAの Φp (Q) に相当する量であることがわかる．
実際の計算は各原子の位置座標に原子質量の平方根をかけた，質量重み付けした座標を用いて行った．上位主成分の数 Nc は 10 とした．これらの主成分の全分散への寄与は 38%であった．t1 = 0.8 ps, t2 = 30.0 ps, τ = 15.0 ps と選び PCRMA を適用した．得られた緩和率は小さいの分布を示している．p = 1 − 10 の緩和時間
順に番号をつける．図 8 は得られた緩和時間 λ−1pは p = 11 − 20 の緩和時間と比べて約 1 桁長い．2 つの発展時間を用いた PCRMA を適用したのは，速い緩和モードの寄与を減らすためであり，p = 11 − 20 の緩和モードはそのような速い緩和モードで，ここでは重要ではない．
（5.16）式を用いて，主成分の自己相関関数 Φp (t)Φp (0) を計算し，シミュレーションで直接計算した結果と比較した．図 9 に p = 1, 2 の場合の結果を示す．ここで示したのは短時間領域

統計数理

236

第 62 巻

第 2 号 2014

図 9．主成分の自己相関関数 Φp (t)Φp (0) の t 依存性．p = 1 の場合を（a）に，p = 2 の場合（5.16）式を用を（b）に示した．実線がシミュレーションから直接計算した結果，破線がいて PCRMA から計算した結果を表す．Reprinted with permission from T. Nagai,
A. Mitsutake, and H. Takano, J. Phys. Soc. Jpn. 82, 023803（2013）. Copyright
2013, The Physical Society of Japan.

であるが，PCRMA の結果がシミュレーションによる直接計算の結果と良く一致していることがわかる．同様の良い一致は p = 3 − 10 でも見られた．このことから，PCRMA の結果は妥当なものと考えられる．
（a）
図 10 に第 1，第 2，第 3 モードベクトルを PCA の場合と PCRMA の場合に示した．図 10
でわかるように，PCA の第 1 モードは外に出ているターンの部分の Pro103 と Gln104 の残基の揺らぎが局所的に大きい．一方，RMA の第 1 モードにはそのような大きな揺らぎは Pro103と Gln104 の近傍には無く，協同的なモードになっている．Pro103 と Gln104 の大きな揺らぎは
RMA の第 3 モードに見ることができる．
ここでの結果が例示するように，PCA は速くて局所的な大きな揺らぎを上位の主成分として取り出すことがある．PCA と PCRMA のモードベクトルの類似性を調べてみると，揺らぎの大きさが必ずしも運動の遅さに対応しておらず，PCRMA が遅い緩和を取り出すのに有効であることがわかる．
6.

まとめ

詳細釣り合いの成り立つシミュレーション法に対して，時間発展演算子の固有値，固有関数から，扱っている系の緩和率，緩和モードをそれぞれ定義することができることを説明した．計算機シミュレーションの時系列データから緩和率分布と緩和モードを変分問題に基づいて評価する方法である緩和モード解析の方法について説明し，孤立線状高分子系の分子動力学シミュ
レーションへの応用例を示した．生体高分子系への応用として，真空中のペプチド系のモンテ
カルロ・シミュレーションと水中のタンパク質系の分子動力学シミュレーションへ適用した例を紹介した．ペプチド系の例では，主成分分析では抽出できない遅い運動を抽出できることを示した．また，緩和モード解析による自由エネルギー地形では，極小点が軸に沿って並ぶ傾向があり，安定・準安定状態間の遷移を良く記述できていると考えられることを見た．タンパク質系の例では，全原子の大きな自由度を系統的に逓減する方法として，主成分分析を先に行い，上位の主成分だけを用いて緩和モード解析を行う，主成分緩和モード解析の方法を紹介した．
さらに，1 つの発展時間だけでは残ってしまう速い緩和モードの影響を除くための方法として，
2 つの発展時間を用いる緩和モード解析の方法を紹介した．2 つの発展時間を用いる主成分緩和
モード解析により，全原子を扱いながら遅い緩和モードを系統的に抽出できることを示した．

生体高分子系の緩和モード解析

237

1/2

図 10．平均構造上の Cα , Cβ , Cγ 原子の位置でのモードベクトル．PCA で得られた Λ1 F 1
1/2
1/2
を
（a）に，Λ2 F 2 を（c）に，Λ3 F 3 を（e）に，PCRMA で得られた h1 を（b）に，h2を（d）に，h3 を（f）に示した．質量重み付けを外し，振幅を 8 倍している．PCA のを数値とし結果には分散 Λn を，PCRMA の結果には振幅 |hp |2 と緩和時間 λ−1pて示してある．Pro103 と Gln104 の Cα 原子を黄色い球で示した．Reprinted withpermission from T. Nagai, A. Mitsutake, and H. Takano, J. Phys. Soc. Jpn. 82,
023803（2013）. Copyright 2013, The Physical Society of Japan.

以上，生体高分子系のシミュレーションから得られる時系列データに対して，緩和モード解析の方法が適用できることを示した．近年，タンパク質の分子動力学シミュレーションが盛んに行われてきているため，その時系列データの動的解析方法として，2 つの発展時間を用いる
主成分緩和モード解析は特に有用であると考えられる．
付A.

録

固有値方程式と変分問題（3.3）式に対し，δR = 0 が固有値方程式（2.13）式，（2.14）式を与えることを見る．δR は次の

統計数理

238

第 62 巻

第 2 号 2014

ように計算される．
(A.1)

⎧
⎫
⎡
⎨
⎬

1
⎣
δR =
δφn (Q)Tτ (Q|Q )ψn (Q ) − Rψn (Q)⎩ ⎭Q φn (Q)ψn (Q)
Q
Q
⎫
⎧
⎤
⎬
⎨
+
φn (Q)Tτ (Q|Q ) − Rφn (Q ) δψn (Q )⎦ .
⎭
⎩

Q

Q

ここで，ψn (Q) = φ̂n (Q)Peq (Q) = φn (Q)Peq (Q) であるから，上式は
⎧
⎫
⎡
⎨
⎬

1



⎣
(A.2) δR =
δφn (Q)Tτ (Q|Q )φn (Q )Peq (Q ) − Rφn (Q)Peq (Q)⎩ ⎭
Q φn (Q)ψn (Q)
Q
Q
⎫
⎧
⎤
⎬
⎨
φn (Q)Tτ (Q|Q ) − Rφn (Q ) δφn (Q )Peq (Q )⎦
+
⎭
⎩
Q
Q
⎧
⎫
⎡
⎨
⎬

1



⎣
=
δφn (Q)Tτ (Q|Q )φn (Q )Peq (Q ) − Rφn (Q)Peq (Q)⎩ ⎭
Q φn (Q)ψn (Q)
Q
Q
⎫
⎧
⎤
⎬
⎨
φn (Q)Tτ (Q |Q)Peq (Q) − Rφn (Q )Peq (Q ) δφn (Q )⎦
+
⎭
⎩
Q
Q
⎧
⎫
⎡
⎨
⎬

1
⎣
=
δφn (Q)Tτ (Q|Q )φn (Q )Peq (Q ) − Rφn (Q)Peq (Q)⎩ ⎭
Q φn (Q)ψn (Q)
Q
Q
⎧
⎫⎤
⎨
⎬

+
δφn (Q )Tτ (Q |Q)φn (Q)Peq (Q) − Rφn (Q )Peq (Q ) ⎦
⎩
⎭
Q
Q
⎧
⎫
⎨
⎬

2
δφn (Q)Tτ (Q|Q )φn (Q )Peq (Q ) − Rφn (Q)Peq (Q)=⎩ ⎭
Q φn (Q)ψn (Q)Q

Qと変形され，δR = 0 より固有値方程式（2.14）式Tτ (Q|Q )φn (Q )Peq (Q ) − Rφn (Q)Peq (Q) = 0
Qが得られ，このときの R の値が固有値 exp(−λn ) となる．ここで，（3.16）式と（2.12）式を用い
式はた．同様に，
（A.2）
(A.3)

⎧
⎫
⎡
⎨
⎬

1


⎣
δR =
δφn (Q)Tτ (Q |Q)φn (Q )Peq (Q) − Rφn (Q)Peq (Q)⎩ ⎭
Q φn (Q)ψn (Q)
Q
Q
⎫
⎧
⎤
⎬
⎨
φn (Q)Tτ (Q|Q ) − Rφn (Q ) δφn (Q )Peq (Q )⎦
+
⎭
⎩
Q
Q
⎫
⎧
⎡
⎬
 ⎨
1


⎣
φn (Q )Tτ (Q |Q) − Rφn (Q) δφn (Q)Peq (Q)=⎭⎩ 
Q φn (Q)ψn (Q)QQ生体高分子系の緩和モード解析

239

⎫
⎤
⎬
+
φn (Q)Tτ (Q|Q ) − Rφn (Q ) δφn (Q )Peq (Q )⎦
⎭
⎩
Q
Q
⎫
⎧
⎬
⎨
2
=
φn (Q)Tτ (Q|Q ) − Rφn (Q ) δφn (Q )Peq (Q )⎭⎩Q φn (Q)ψn (Q)

Q
⎧
⎨

Q

と変形され，δR = 0 より固有値方程式（2.13）式
φn (Q)Tτ (Q|Q ) − Rφn (Q ) = 0
Qが得られ，このときの R の値が固有値 exp(−λn ) となる．
φn (Q) = φ̂n (Q) となるモンテカルロ法や Brown 動力学法では，緩和率 λn は全て 0 以上の実数で，δR = 0 となる R の値は極値となり，変分問題 δR = 0 は通常の変分法としての意味をもつ．これに対し，分子動力学法では，緩和率 λn が複素数となるため，δR = 0 となる R の値は一般に極値とはならず，通常の変分法としての意味はもたない．1 次の変分 δR = 0 が固有値方
程式を与え，そのときの R の値が固有値になるという意味で，R は有用な量である．
Hirao et al.（1997）
では，
（A.1）
式の δφn (Q) と δψn (Q) があたかも独立であるかのように扱っている．得られる結果は正しいが，詳細釣り合いの条件（2.11）式が成立しているので，δφn (Q)と δψn (Q) は独立ではなく，上記のように扱わなくてはならない．
B.

試行関数の選び方

式の選び方について考える．緩和モード解析により，複数の自由度がどのよう試行関数（3.5）に協力して遅い緩和を実現しているかを明らかにすることができる．従って，対象としている
運動を構成している自由度のうち，重要と考えられる自由度の線形結合を試行関数とする．例
えば，ランダムスピン系でスピンの協力的な運動を取り出すには，各スピンの値の線形結合を．第 3 節の孤立高分子系の例では，高分子の重心周りでの用いる（Takano and Miyashita, 1995）
変形運動を問題にしたので，各セグメントの重心からの相対座標の線形結合を用いた．第 4 節，
第 5 節で出てくる生体高分子系の例では，並進と回転を除いた後の平均的な構造のまわりの遅い揺らぎを調べるため，平均的な構造からの原子の変位の線形結合を用いている．遅い緩和に重要だと思われる自由度が予想できていれば，それらの線形結合を用いれば，より良い試行関数を構成することができる．
次に，試行関数の中の t0 の役割について考える（Mitsutake et al., 2011）．簡単のため，純粋（4.1）
式を用いる．
（3.5）式に相当する試行関数（近似的緩緩和だけの場合を考え，第 4 節の変数
和モード）は(B.1)
Xp (Q) =fp,i Ri (t0 /2; Q)iとなる．Ri (Q) が厳密な緩和モード {φn (Q)} で(B.2)
Ri (Q) =ri,n φn (Q)nと展開されるとすると，
（3.6）
式より，


(B.3)
Ri (t; Q) =ri,n φn (Q )Tt (Q |Q) =ri,n exp(−λn t)φn (Q)
Q

n

n

統計数理

240

第 62 巻

第 2 号 2014となる．ここで，λn は φn (Q) の緩和率で，番号は λn の小さい順につけるものとする．第 4 節の例では，R が 3N − 6 個の自由度をもち，一般化固有値問題（4.8）式および（4.9）式を解くこと
で，R に含まれる遅い方から 3N − 6 個の緩和モードを近似的に推定する．もし，R が 3N − 6個の緩和モードしか含んでいなければ，t0 = 0 としても緩和モードを精度良く推定することができる．しかし，R には最も遅い 3N − 6 個の緩和モード以外に，より速い緩和モードが多数
式を見ると，速い緩和モードの遅い緩和モードに対する相対的寄与の大き含まれている．（B.3）
さ，例えば，1 ≤ m ≤ 3N − 6 < n として，n 番目の緩和モードの m 番目の緩和モードに対する相対的寄与の大きさは |ri,n /rj,m | exp[−(λn − λm )t] で与えられ，λn > λm であることから，tを大きくすることで，速い緩和モードの相対的寄与が小さくなることがわかる．（B.1）式の中のt0 は，試行関数に含まれる推定範囲外の速い緩和モードの寄与を減らすことで，緩和モード推定の精度を上げる役割をしている．
謝辞

共同研究者である宮下精二氏，平尾秀友氏，古関幸子氏，光武亜代理氏，飯島寛允氏，長井俊樹氏，山本卓氏，小泉祐太氏，そして，この報告を書く機会を与えて下さいました伊庭幸人
氏に深く感謝致します．

参

考

文

献

Doi, M. and Edwards, S. F.（1986）. The Theory of Polymer Dynamics, Oxford University Press, Oxford.
Drozdov, A. N.（2000）. Variational theory of activated rate processes for an arbitrary barrier, PhysicalReview E, 62, 1879–1892.
Hirao, H., Koseki, S. and Takano, H.（1997）. Molecular dynamics study of relaxation modes of a single
polymer chain, Journal of the Physical Society of Japan, 66, 3399–3405.
Koseki, S., Hirao, H. and Takano, H.（1997）. Monte Carlo study of relaxation modes of a single polymer
chain, Journal of the Physical Society of Japan, 66, 1631–1637.
Kremer, K. and Grest, G. S.（1990）. Dynamics of entangled linear polymer melts: A molecular-dynamics
simulation, The Journal of Chemical Physics, 92, 5057–5086.
Mitsutake, A., Hansmann, U. H. E. and Okamoto, Y.（1998）. Temperature dependence of distributions of conformations of a small peptide, Journal of Molecular Graphics and Modelling, 16,
226–238.
Mitsutake, A., Iijima, H. and Takano, H.（2011）. Relaxation mode analysis of a peptide system: Comparison with principal component analysis, The Journal of Chemical Physics, 135, 164102-1-15.
Nagai, T., Mitsutake, A. and Takano, H.（2013）. Principal component relaxation mode analysis of an
all-atom molecular dynamics simulation of human lysozyme, Journal of the Physical Society ofJapan, 82, 023803-1-4.
Risken, H.（1996）. The Fokker-Planck Equation: Methods of Solution and Applications, SpringerVerlag, Berlin, Heidelberg.
Takano, H. and Miyashita, S.（1995）. Relaxation modes in random spin systems, Journal of the PhysicalSociety of Japan, 64, 3688–3698.
Talkner, P.（1994）. Finite barrier corrections for the Kramers rate problem in the spatial diﬀusion
regime, Chemical Physics, 180, 199–213.

Proceedings of the Institute of Statistical Mathematics Vol. 62, No. 2, 221–241 (2014)

241

Relaxation Mode Analysis of Biopolymer Systems
Hiroshi Takano
Faculty of Science and Technology, Keio University

For simulation methods that satisfy the detailed balance condition, relaxation rates
and modes of a system concerned can be defined from the eigenvalues and eigenvectorsof the time-evolution operator, respectively. The relaxation mode analysis method, which
estimates the relaxation rate distribution and relaxation modes on the basis of the variational method, is explained. Applications of the method to Monte Carlo simulation of a
peptide system and molecular dynamics simulation of a protein system are reviewed.

Key words: Relaxation modes, relaxation rates, relaxation mode analysis, simulations, time-displacedcorrelation functions.