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

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


生体分子シミュレーションの摂動解析
小山 洋平†（受付 2013 年 12 月 20 日；改訂 2014 年 7 月 18 日；採択 7 月 22 日）

要旨タンパク質の機能を原子レベルから詳細に理解するためにはタンパク質主鎖や側鎖の異なる構造状態の同定，化合物・DNA・他のタンパク質との様々な結合状態の同定，これらに対する
水やイオンなどの溶媒の影響，について総合的に評価する必要がある．このような複雑な構造
状態とその相互作用を明らかにするために，原子間相互作用に対する摂動解析を用いた我々の取り組みについて解説する．まず，タンパク質内あるいは複合体内相互作用に対する摂動解析である PEPCA（potential energy principal component analysis）を導入する．PEPCA は相互作
用エネルギーに対する主成分分析を実行することで，その主成分得点によりタンパク質内あるいは複合体内の相互作用パターンが異なる構造状態が分離され，固有ベクトル成分から構造状
態間の差異を生み出している相互作用の情報が得られる．次に，タンパク質または複合体と溶媒間相互作用の摂動解析である DIPA（distance-dependent intermolecular perturbation analysis）
を導入する．DIPA は原子間力と平均溶媒原子数の積に対して関数主成分分析を実行することで，その主成分得点により溶媒との相互作用パターンが異なる構造状態が分離され，固有関数
成分により構造状態間の相違を生み出している溶媒との距離依存的な相互作用の情報が得られる．最後に，10 個のアミノ酸からなるシニョリンのフォールディングシミュレーションに対して PEPCA と DIPA を適用し，固有ベクトル成分または固有関数成分と主成分得点の散布図を同時に表示する可視化法であるバイプロットを用いた結果の解釈について解説する．
キーワード：分子動力学シミュレーション，摂動解析，主成分分析，関数主成分分析，
バイプロット．
1.

はじめに

生体中には様々な種類の分子が存在しているが，タンパク質が触媒反応などの機能を発現するためには特定の分子を認識し，その機能を選択的に発現する必要があると考えられる（特異性）．
このことからタンパク質に様々な刺激を加え，高い応答性を示す刺激を調べることでタンパク質の機能発現メカニズムや特異性について明らかにできる可能性がある．Ming らはタンパク質構（normal mode anal造と化合物構造の結合前後での同時分布（joint distribution）を基準振動解析
ysis）で解析的に求めた．化合物の結合に伴うタンパク質構造の応答性を定量化するために，結合前後での同時分布変化の大きさをカルバック・ライブラーダイバージェンス（Kullback-Leibler
を用いて測った．化合物をタンパク質表面の様々
divergence; 相対エントロピー; Kullback, 1959）な部位に発生させ，各部位でカルバック・ライブラーダイバージェンスを計算すると，実際の化合物の結合部位で大きな分布変化が生じることを見出した（Ming and Wall, 2005）．その後，†理化学研究所 生命システム研究センター：〒 650–0047 兵庫県神戸市中央区港島南町 1–6–5

186

統計数理

第 62 巻

第 2 号 2014

同時分布ではなく化合物発生前後でのタンパク質の構造分布（周辺分布）
を基準振動解析により求め，タンパク質の構造分布変化の大きさをカルバック・ライブラーダイバージェンスにより定量化した．様々なタンパク質と化合物のデータセットを用いて，実際の化合物の結合部位ではタンパク質構造の分布変化が大きくなり，化合物結合部位の予測に有用であることを示した．
（Ming and Wall, 2006）
通常の機械とは異なり，タンパク質のような生体分子機械は自発的な動きや周りの分子との相互作用を通じて常にゆらいでいる．分子動力学（MD）シミュレーションを用いてこのような生体分子や溶媒のゆらぎを詳しく調べることが可能である．平衡統計力学では平衡状態に摂動が加わった新たな平衡状態でのゆらぎを，摂動前の平衡状態に摂動法を適用することで推測することが可能である．特に，摂動が “微小” であり，摂動ハミルトニアンが時間にあらわに依存しないときには，摂動前後での物理量の平均値の差が未摂動平衡状態での物理量と摂動ハミルトニアンとの共分散で表されるという静的な線形応答理論（static linear response theory）が成．池口らは化合物の結合を原子に一定力を加えるという摂動でモデル化立する（Zwanzig, 2001）
し，化合物が非結合状態での MD シミュレーションに対して静的な線形応答理論を適用することで，この摂動による平均構造変化の予測を行った（Ikeguchi et al., 2005）．この結果，実際の化合物結合前後での構造変化と良く一致すること，力の向きを多少変化させても構造変化の方向はあまり変わらないという頑強性があることを示した．
タンパク質などの生体分子の機能はその構造形成（フォールディング）や（部分的な）アンフォー
ルディング，化合物・DNA・その他のタンパク質との様々な部位での結合・解離，それに前後する構造変化などが周囲の水やイオンなどの溶媒の影響を受けながら生じる．従って，生体分子の機能を原子レベルから詳細に理解するためにはこのような複雑な構造状態の変化を体系的に取り扱う必要がある．生体分子構造の複雑性は生体分子内および生体分子・溶媒分子間におけるファンデルワールス相互作用や静電相互作用などの原子間相互作用により生じており，複数の構造状態が存在する場合には，これらのパターンが異なっていると考えられる．このため，
本稿ではどのような構造状態が存在するかを原子間相互作用により特徴付け，各状態間の差異を生み出している影響の大きな相互作用を明らかにし，それらの各状態への寄与を明らかにすることを目指す．このためには，相互作用の生体分子構造状態への影響を総合的に評価する必要がある．生体分子構造の静的なゆらぎの情報は生体分子の構造分布から得ることができるた
め，相互作用の強弱を変化させたときに生体分子の構造分布がどれほど変化するかを調べることでその影響を評価することにする．構造状態間の差異に関係していない相互作用の強弱を変化させても生体分子の構造分布はあまり大きく変化しないと考えられる一方で，差異を生み出している相互作用の強弱を変化させたときには大きな構造分布変化が生じると予想される．従って，相互作用の強弱を変化させ，大きな構造分布変化を生み出す相互作用の組み合わせを探索することで，生体分子の構造状態に大きな影響を及ぼす相互作用を見出すことができると考えられる．また，これらの相互作用の各状態への影響は，相互作用の強弱を変化させたときの構造状態の比率変化により特徴付けられると予想される．
相互作用の強弱を変化させることは摂動ポテンシャルエネルギーを導入することで実現でき，
このような摂動により生体分子の構造状態を理解しようとする上記のアプローチを本稿では相互作用の摂動解析と呼ぶことにする．摂動に伴う生体分子構造の分布変化はカルバック・ライ
ブラーダイバージェンスを用いて定量化を行う．フォールディング，部分的なアンフォールディ
ング，複合体の形成・解離，といったような基準振動解析が適用できないような大きな構造変化も扱えるように，平衡状態の MD シミュレーションの結果に対して摂動法を適用することでカ
ルバック・ライブラーダイバージェンスを推定する．生体分子内相互作用の摂動解析（Koyamaと生体分子・溶媒分子間相互作用の摂動解析（Koyama et al., 2011）について紹介を
et al., 2008）

生体分子シミュレーションの摂動解析187行った後，これらの摂動解析の結果がバイプロットを用いて理解できることを示す．
2.

摂動解析

2.1 摂動による生体分子構造の分布変化
生体分子の座標を q ，溶媒分子の座標を q  ，系のポテンシャルエネルギーを V (q, q  ) とする.
ここで，生体分子は複数種類あってもよい．原子数・体積・温度一定（N V T ）条件での平衡状態における分布は逆温度を β とするカノニカル分布 ρ(q, q  ) = e−βV (q,q ) /Z で表される．Z は規格化定数に相当する分配関数である．摂動ポテンシャルエネルギー ΔV (q, q  ) を用いて系に摂動を与えることを考える．摂動後の系のポテンシャルエネルギーは V  (q, q  ) = V (q, q  ) + ΔV (q, q  )となり，摂動後の分配関数を Z  とすると，摂動後の分布は


(2.1)

ρ (q, q  ) =1 −βV  (q,q  )e−βΔV (q,q )e=
ρ(q, q  )Z
e−βΔV となる. A は未摂動平衡状態での物理量 A(q, q  ) の平均を表す．このとき，摂動後の生体分子の構造分布は e−βΔV  q




(2.2)
ρ (q) = ρ (q, q )dq =ρ(q)e−βΔV と書ける．ここで，ρ(q) は摂動前の生体分子の構造分布，A|q は生体分子構造 q を固定した状
態での物理量 A(q, q  ) の未摂動平衡状態での条件付期待値を表す．摂動による分布変化をカル
バック・ライブラーダイバージェンス D(ρ (q)||ρ(q)) ≡ ρ (q) ln ρ (q)/ρ(q)dq で定量化し，式（2.2）を代入すると
1(2.3)
D(ρ (q)||ρ(q)) = var(βΔV |q) + · · ·2と展開できることが分かる．var(A) は未摂動平衡状態における物理量 A(q, q  ) の分散を表す．
以上のようにして，摂動ポテンシャルエネルギー ΔV による生体分子構造の分布変化は ΔV の条件付期待値の分散で近似できることが分かる．
ΔV が生体分子構造 q だけに依存する場合には ΔV |q = ΔV (q) となり，条件付期待値を計算する必要が無くなる．生体分子内の相互作用に摂動を加える場合がこれに該当する．一方で，生体分子と溶媒分子の間の相互作用に摂動を加える場合には ΔV が q と q  に依存するため，条件付期待値を具体的に評価する必要が生じる．このように，生体分子内相互作用への摂動と生体分子・溶媒分子間相互作用への摂動では，式（2.3）の条件付期待値に対する扱いが変わるため，
以下ではこの二つを分けて考える．
2.2 生体分子内相互作用の摂動解析（PEPCA）
生体分子構造の複雑性はファンデルワールス相互作用や静電相互作用などの原子間相互作用の組み合わせにより生じている．生体分子構造への影響が大きな生体分子内相互作用を明らかにするために，その強弱を変化させる摂動解析を考える．ある生体分子内相互作用のポテンシャル
エネルギーを Vk (q) とするとき，摂動パラメーター λk を導入してエネルギー項を (1 + λk )Vk (q)とすれば，λk > 0 で相互作用を強め，λk < 0 で弱めることができる（図 1（a））．これより，相互作用の強弱を変化させるには λk Vk (q) という摂動項を加えれば良いことが分かる．生体分子
内の M 個の相互作用 V (q) = (V1 (q), . . . , VM (q)) の強弱を変える摂動は M 個の摂動パラメー
ター λ = (λ1 , . . . , λM )T を用いて摂動ポテンシャルエネルギーを

188

統計数理

第 62 巻

第 2 号 2014

図 1．生体分子内相互作用の摂動解析（PEPCA）．

図 2．等価原子の番号置換に対するポテンシャルエネルギーの不変性と摂動による影響.

(2.4)

ΔV (q) =M

λi Vi (q)i=1とすることで実現できる．
アスパラギン酸の側鎖に含まれるカルボキシル基 CO2 の酸素原子やアルギニンのグアニジ（図 2（a））に対して，ノ基に含まれる NH2 の水素原子のように，等価な原子の番号の入れ替え（b）; Malolepsza et al.,
ポテンシャルエネルギーが不変であるという置換対称性が存在する（図 2．このような等価な原子との相互作用に対して異なる摂動パラメーターを用いて摂動を加2010）（c））．一方で，えると，摂動後のポテンシャルエネルギーの置換対称性が失われてしまう（図 2
等価な原子との相互作用のどちらか一方が生体分子の構造分布により大きな影響を及ぼすと決めることはできないため，これらの相互作用に対する摂動パラメーターには同一のパラメーターを用いることにする．これは置換対称性のあるポテンシャルエネルギーのグループをまとめて
Vk と再定義することに相当し，このとき，摂動後もポテンシャルエネルギーの置換対称性が保（d））．
たれることが分かる（図 2という条件で生体分子内相互作用の強弱を様々に変化させ摂動の強さが一定（|λ| = δ, δ > 0）
たときに，大きな構造分布変化を示す相互作用の摂動を調べることで，生体分子構造に対して

生体分子シミュレーションの摂動解析

189

大きな影響を及ぼす生体分子内相互作用が明らかにできると考えられる（図 1（b））．生体分子内による生体分子構造の分布変化の大きさ（2.3）は相互作用の摂動
（2.4）

M

1
(2.5)
D(ρ (q)||ρ(q)) = varλi βVi (q) + · · ·
2
i=1
(2.6)

=

1 T
λ cov(βV )λ + · · ·2と書ける．ここで，cov(βV ) は未摂動平衡状態での βV (q) の共分散行列を表す．|λ| = δ というを最も大きく変化させるような λ を探すためには，δ が小さく，式（2.5）
条件で構造分布変化（2.5）の分散を最大化すればよいことが分かる．式（2.5）の近似が成り立つという条件において式（2.5）の分散内を係数 λi と統計量 βVi (q) の線形結合とみなすと，|λ| = δ という条件での分散の最（principal component analysis; PCA; Jolliﬀe,
大化は統計量 βV (q) または −βV (q) の主成分分析
により実行できる．従って，共分散行列 cov(βV ) の固有値を降順に σi2 ，対応する固有ベ
2002）
クトルを ui とすると，λ = δu1 のときに |λ| = δ という条件で分散が最大化されることが分かる．同様にして，λ = δui は |λ| = δ および uT1 λ = 0, . . . , uTi−1 λ = 0 という条件での分散の最大化に対応する．
摂動 λ での生体分子構造の分布を ρλ (q) とする．摂動 λ = δui による生体分子構造の分布変は化の大きさ（2.6）
1(2.7)
D(ρδu i (q)||ρ(q)) ≈ δ 2 σi22となり，共分散行列の固有値が分布変化の大きさに対応することが分かる．これより，大きい
固有値に対応する固有ベクトル（上位の固有ベクトル）
により生体分子構造に大きな影響を及ぼす生体分子内相互作用が得られると考えられる．また，gi (q) ≡ −βuTi (V (q) − V ) と定義する
（2.2）
（2.4）
，
を用いて
と，構造 q の比率変化は式
(2.8)

ln

ρδu i (q)= δgi (q) + · · ·ρ(q)と書ける．ここで導入した gi (q) はその定義から −βV (q) を用いた PCA の第 i 主成分得点であにより主成分得点 gi (q) の符号で摂動による構造比率の増減を近似的ることが分かる．式（2.8）に推定でき，対応する固有ベクトルの相互作用がどのような構造状態に対して影響を及ぼしているかを推測することができる．この影響を可視化する手法として次節でバイプロットを用いて解説する．以上から，ポテンシャルエネルギー −βV (q) の PCA（potential energy PCA; PEPCA;を行うことで生体分子内相互作用の摂動解析が実行
Koyama et al., 2008; Koyama et al., 2011）（c））．
できることが分かる（図 1
2.3 バイプロット
ここでは固有ベクトル成分の散布図と主成分得点の散布図を同時に表示する可視化法である
バイプロット（biplot; Gabriel, 1971; Gower and Hand, 1996）を用いることで PEPCA 固有ベク
トルやその線形結合による摂動が引き起こす構造比率の変化を調べることが可能であることについて説明する．以下では 2 次元のバイプロットで説明を行うが 3 次元のバイプロットでも同様の解釈ができる．PEPCA 固有ベクトル ui の k 番目の成分を Uki とする．第 i, j 固有ベクトと第 i, j 主成分得点の散布図（gi (q), gj (q)）を同時に表示した 2 次元バ
ル成分の散布図（Uki , Ukj ）．
イプロットを考える
（図 3）
（図 3
（a））．摂動 λ = δui により k
最初に第 i 固有ベクトルの摂動による影響について考える
番目の相互作用のポテンシャルエネルギーは (1 + δUki )Vk (q) となる．Uki > 0 であれば相互作

190

統計数理

第 62 巻

第 2 号 2014

図 3．バイプロット. 摂動による構造比率の増加・減少は式（2.8）の近似に基づく．

用が強まり，Uki < 0 であれば弱まる．また，式（2.8）の近似の下で gi (q) > 0 であれば構造比率が増加し，gi (q) < 0 であれば減少する．このようにして，バイプロット上で左側の相互作用を弱め，右側の相互作用を強めると左側の主成分得点に対応する構造の比率が減少し，右側の構（a1））．同様にして，摂動 λ = −δui を考えることで左造の比率が増加することが分かる（図 3側の相互作用を強め，右側の相互作用を弱めると，左側の構造の比率が増加し，右側の構造の（a2））．第 j 固有ベクトルの摂動による構造比率変化につい比率が減少することが分かる
（図 3
（b））．
ても同様に，上記の左右を下上に置き換えればよい（図 3
また，固有ベクトルの線形結合 ui,j (θ) ≡ cos θui +sin θuj , −π ≤ θ < π による摂動 λ = δui,j (θ)
を導入することで左右や下上で分離される領域の構造比率変化だけではなく，原点を通り単位法
線ベクトル n ≡ (cos θ, sin θ) を持つ直線により分離される領域についても同様の考え方ができ（c））．ui,j (θ) の k 番目の成分は n · (Uki , Ukj ) と書け，固有ることが次のようにして分かる（図 3を軸 n へ射影することで得られる．また，摂動 λ = δui,j (θ) による構造
ベクトル成分（Uki , Ukj ）（，2.4）
を用いると g(q) ≡ cos θgi (q) + sin θgj (q) と定義することにより式（2.8)
比率変化は式（2.2）と同様の式が成り立つことが分かる．g(q) = n · (gi (q), gj (q)) と書けることから g(q) は主成分得
点（gi (q), gj (q)）を n に射影することで得られる．以上から，固有ベクトル成分（Uki , Ukj ）を軸 nに射影した値が正である相互作用を強め，負である相互作用を弱めると主成分得点（gi (q), gj (q)）
を軸 n に射影した値が正である構造の比率が増加し，負である構造の比率が減少することが分（c1））．
かる（図 3以上の結果を利用することにより，バイプロット上で構造状態間の差異を生み出している相（a）において相互作用 k = 5, 6, 7, 8 を強め，互作用を明らかにすることができる．例えば，図 3
相互作用 k = 1, 2, 3, 4 を弱めると状態 C の比率が増加し，状態 A の比率が減少することが分かる．このことから相互作用 k = 5, 6, 7, 8 は状態 C を状態 A に対して相対的に安定化（あるいはしている相互作用であると推測できる．なお，上状態 A を状態 C に対して相対的に不安定化）
記のようにバイプロット上で左側の相互作用を強めると左側の構造比率が増加するというよう生体分子シミュレーションの摂動解析

191

図 4．等価分子の番号置換に対するポテンシャルエネルギーの不変性と摂動による影響.

図 5．距離に依存した生体分子・溶媒分子間相互作用の摂動解析（DIPA）．

な解釈をするためには βV (q) ではなく −βV (q) を用いて PCA を行うことが重要である．
2.4 距離に依存した生体分子・溶媒分子間相互作用の摂動解析（DIPA）
次に，水やイオンなどの溶媒分子が生体分子に及ぼす影響を明らかにするために，生体分子・
溶媒分子間相互作用への摂動解析を考える．まず，2.2 節で導入した相互作用の強弱を変化させる摂動を適用した場合について考える．溶媒分子の場合には分子内の置換対称性だけではな（a））
に対するポテンシャルエネルギーの置換対称性が存在す
く，溶媒分子の番号入れ替え
（図 4
（b）
）
．このため，摂動後の置換対称性を保つためには置換対称性のある相互作用を全てる（図 4の溶媒分子についてまとめて一つのポテンシャルエネルギー項として扱う必要があるが，このポテンシャルエネルギー項に対して摂動を加えると，生体分子近くにある溶媒分子も遠くにある溶媒分子も同じように相互作用の強さが変化する．このように，溶媒分子の距離によらず一様に相互作用の強さを変化させる摂動解析（intermolecular perturbation analysis; IPA; Koyama
et al., 2011）
では主成分得点による構造状態の分離が不十分であり，溶媒の生体分子に対する距
離依存的な寄与を考慮する必要があることが示唆された．
溶媒分子の距離依存的な寄与を考えるためには，距離に応じて相互作用の強弱を変えることがで
（a）
（
，b）
）
．原子間距離 r での相互作用のポテンシャルエネルギー
きる摂動解析を行えばよい（図 5（1 + λ(r)）倍となるを φ(r)，距離に依存した摂動パラメーターを λ(r) とすると原子間に働く力が

192

統計数理

第 62 巻

第 2 号 2014

ような摂動を考えたい．このような摂動に対応するポテンシャルエネルギーは F (r) ≡ −dφ(r)/drとすると⎧⎨ rc λ(r )F (r )dr , r ≤ rcr
(2.9)
φ(r; λ) ≡⎩ 0, r > rcにより実現できる．ここで，カットオフ距離 rc を導入した．rc の値を大きくしたときに摂動解
析の結果が収束するかどうかを確認し，収束していればカットオフが無い（rc = ∞）場合の結果を表していると考えられる．また，摂動パラメーターが距離に依存しない場合と異なり，摂動項が λ(r)φ(r) とはならないことに注意が必要である．距離に依存した摂動の場合にも個別に摂（c）），置換対称性のある相互作用につい動を加えると摂動後の置換対称性が失われるため（図 4（d））．
ては同一の摂動パラメーターを用いて摂動を加える必要がある（図 4以上の結果に基づき，置換対称性を持つ相互作用のグループ k に対して，置換対称性を保つ
距離依存的な摂動は，共通の相互作用のポテンシャルエネルギー関数を φk (r)，共通の摂動パラメーターを λk (r) とすると
 
(2.10)ΔVk (q, q  ) ≡φk (rij ; λk )i∈Ik j∈Jk

で与えられる．Ik は相互作用 k に属する生体分子の等価原子の集合，Jk は相互作用 k に属する
溶媒の等価原子の集合，rij は原子 i と原子 j の間の距離を表す．次に，摂動による生体分子構にある条件付期待値 ΔVk |q の評価について考える．生体分子の構造 q が固造の分布変化（2.3）
定されているとき，生体分子の原子 i から距離 r にある Jk に含まれる溶媒原子の平均数密度を
niJk (r|q) とすると式（2.10）の条件付期待値は  rc(2.11)
ΔVk |q =φ(r; λk )niJk (r|q)dri∈Ik0と書ける．ここで，生体分子の構造 q が固定されているときの生体分子の原子 i から距離 r 以
r
内にある平均溶媒原子数 NiJk (r|q) ≡ 0 niJk (r |q)dr および(2.12)
fk (r|q) ≡ Fk (r)NiJk (r|q)
i∈Ikを導入し，式（2.11）の部分積分を行うと(2.13)

ΔVk |q = rc
0

λk (r)fk (r|q)drと書けることが分かる．
以 上 か ら ，生 体 分 子・溶 媒 分 子 間 相 互 作 用 に 対 す る 距 離 に 依 存 し た 摂 動 ΔV (q, q  ) =M（2.3）は，f (r|q) = (f1 (r|q), . . . , fM (r|q))T および
k=1 ΔVk (q, q ) による構造分布変化の大きさ
T
λ(r) = (λ1 (r), . . . , λM (r)) を用いて

 M 
 rc
1
(2.14)D(ρλ (q)||ρ(q)) = var β
λk (r)fk (r|q)dr + · · ·2
k=1 0(2.15)と書ける．ここで，内積 λ, λ  ≡=1
λ, Ĉλ + · · ·
2

 rc
0

λ(r)T λ (r)dr および

生体分子シミュレーションの摂動解析

193

図 6．DIPA 実行手順.

(2.16)Ĉλ(r) ≡ rc0cov(βf (r|q), βf (r |q))λ(r )drのように作用する共分散作用素（covariance operator）
Ĉ を導入した．cov(βf (r|q), βf (r |q)) は
その i, j 要素が cov(βfi (r|q), βfj (r |q)) となるような共分散行列である．これにより生体分子構造の分布を最も大きく変化させるような摂動 λ(r) を見つけるためには式（2.14）の分散を最大化（2.14）の分散の最大化はすれば良い．|λ| ≡ λ, λ と定義すると，|λ| = δ という条件での式（functional PCA; FPCA; Jolliﬀe,
関数データ βf (r|q) または −βf (r|q) を用いた関数主成分分析
により実行できる．従って，共分散作用素 Ĉ の降順の固有
2002; Ramsay and Silverman, 2006）
2
値を σi ，対応する固有関数を ui (r) とすると λ(r) = δu1 (r) のとき |λ| = δ という条件で分散が最大になる．同様にして，λ(r) = δui (r) は |λ| = δ および u1 , λ = 0, . . . , ui−1 , λ = 0 という条件で分散を最大化することが分かる．摂動 λ(r) = δui (r) による構造分布変化の大きさ（2.7）となり，構造比率の変化（2.2）
は gi (q) ≡ −βui (r), f (r|q) − f (r|q) と定義す
（2.15）は式
ると式（2.8）となる．gi (q) は定義から −βf (r|q) を用いた FPCA の第 i 主成分得点である．以上から，−βf (r|q) の FPCA を行うことで距離に依存した生体分子・溶媒分子間相互作用の摂動
解析（distance-dependent intermolecular perturbation analysis; DIPA; Koyama et al., 2011）が実（c））．
行できることが分かる（図 5
DIPA を実行するには平均溶媒原子数 NiJk (r|q) を推定する必要がある．このためには，生体分子を固定して，溶媒分子のみを動かす MD シミュレーションを実行する必要があるが，この（belly）法として知られような特定の自由度だけを動かす MD シミュレーションの方法はベリーている（古明地 他, 2000; Case et al., 2008）．また，生体分子構造 q は平衡状態からサンプリン
グする必要がある．以上から，DIPA の実行手順は次のようになる．まず，生体分子と溶媒分子
（図 6
（a）
）
．この結果得られた生が両方とも動く通常の全原子 MD シミュレーションを実行する
体分子構造を適当な間隔で初期構造として，ベリー法による溶媒分子のみを動かす複数の MD
シミュレーションを実行する（図 6（b））．各シミュレーションで平均溶媒原子数 NiJk (r|q) を推
（図 6
（c）
）
，関数データ行列 −βf (r|q) に対する FPCA
定後，式（2.12）を用いて fk (r|q) を計算しを実行すればよい
（図 6
（d）
）．

194

3.

統計数理

第 62 巻

第 2 号 2014

シニョリンフォールディングシミュレーションの相互作用摂動解析

10 個のアミノ酸（Gly1-Tyr2-Asp3-Pro4-Glu5-Thr6-Gly7-Thr8-Trp9-Gly10）からなり，β ヘアピン構造に折り畳まれるシニョリン（chignolin; Honda et al., 2004）のフォールディングシミュ
レーションを用いて，2.2 節で導入した生体分子内相互作用の摂動解析である PEPCA と 2.4 節で導入した距離に依存した生体分子・溶媒分子間相互作用の摂動解析である DIPA を適用し，バ
イプロットを用いて結果の解釈を行う．
3.1 生体分子内相互作用の摂動解析（PEPCA）の適用（PDB ID : 1UAO）の 1 番目の構造に対して，シニョリンから 20 Å
シニョリンの NMR 構造
以内が水に満たされるように，水分子 6437 個とイオン（12Na+ と 10Cl− ）を追加した．ポテン（1 気圧）
一定（N pT ）シミュレー
シャルエネルギーの極小化を行った後，100 ピコ秒の温度・圧力
ションを行い，1 気圧の密度に体積を緩和させた結果，一辺が 58.4 Å の立方体となった．この後，1 マイクロ秒の温度・体積一定（N V T ）シミュレーションを行った．いずれのシミュレーショ
◦ンもシニョリンの融解温度である 39 C（Honda et al., 2004）で行った．MD シミュレーション
により Langevin 方程式と周期境界条件を用いて実行した．シは AMBER10（Case et al., 2008）
ニョリンは 138 原子からなり，分子内相互作用数は 18, 919 であるが，置換対称性を保つように
グループ化すると M = 12, 927 となった．1 マイクロ秒の構造データから 1 ナノ秒ごとに置換，M × n 中心化データ行列の対称性のあるポテンシャルエネルギーを計算し（n = 1000 データ）を行うことで PCA を実行した（付録 A）．
特異値分解（SVD）
PEPCA の固有値
（図 7）
は第 1, 2 固有値が不連続的に大きく，第 3 固有値以降は連続的になっており，第 1, 2 固有ベクトルでシニョリンの構造状態の安定性を決定している相互作用が特徴（図 8：点）を見ると，大きく分けて付けられることが期待される．第 1, 2 主成分得点の散布図
4 つの状態があることが分かる．このうち，左側の状態に対応する構造は β ヘアピン構造であ，右上の状態に対応する構造は「ほどけた状ることから「正しく折り畳まれた状態」（Native; N）
態」（Unfold; U）と呼ぶことにする．一方，シニョリンのマルチカノニカル MD シミュレーショ
ンの結果により Native 状態である β ヘアピン構造では Asp3H（本稿ではバイプロットでの表が Thr8O と水素結合しているが，Asp3H が Gly7O と水素結合した
記上，3H などと省略する）
Misfold 状態が存在することが見出された（Satoh et al., 2006）．右下の状態に対応する 2 つの構「誤って折り畳まれた状態 1, 2」（Misfold 1, 2; M1, M2）と造はこの Misfold 状態に相当するため呼ぶことにする．

図 7．PEPCA 固有値（Koyama et al., 2011）.

生体分子シミュレーションの摂動解析

195

図 8．PEPCA バイプロット（Koyama et al., 2011）．6OG1-8HG1 などは 6 番目のアミノ（トレオニン）のヒ酸（トレオニン）のヒドロキシル基の酸素原子（O）と 8 番目のアミノ酸
ドロキシル基の水素原子（H）との静電相互作用を表す．小さい点は第 1, 2 主成分得点，小さい四角は第 1, 2 固有ベクトル成分の散布図を表す．第 1 または第 2 固有ベクトル成分の絶対値が上位 20 位までに入る相互作用のみを示す．PEPCA の解析には二面角やファンデルワールス相互作用のポテンシャルエネルギーも含むが，図の相互作用は全て静電相互作用である．また，立体構造において見やすくするために 9 番目のアミノ酸（トリプトファン）の側鎖は表示していない．

2.3 節の議論に基づき，この 4 状態の差異を生み出している相互作用は主成分得点の散布図と固有ベクトル成分の散布図（図 8：四角）を同時に表示したバイプロットを用いて理（図 8：点）解できる．N 状態と M1+M2 状態は左上から右下方向の軸で最もよく分離できる．このため，これらの状態の差異を生み出している相互作用は左上と右下の相互作用を見ればよい．N 状態を安定化あるいは M1+M2 状態を不安定化している相互作用は左上にある相互作用を見ればよく，左上にある 6OG1-8HG1 や 3OD-6HG1 という相互作用は酸素原子 O が負の部分電荷，水素（静電引力）である．従って，これら
原子 H が正の部分電荷を持つため，引力的な静電相互作用（図 8 の vdWの相互作用は N 状態を安定化していることが推測できる．実際に，N 状態の構造
球表示．以下では N 構造などとする．）
で 6OG1 と 8HG1 や 3OD と 6HG1 の間の相互作用が確認できるが，M1 や M2 の構造では見られないことが分かる．また，M1+M2 状態を安定化ある
いは N 状態を不安定化している相互作用は右下の相互作用を見ればよい．右下の 3OD-6OG1 などが集まっている相互作用は全て静電斥力であり，N 状態を不安定化していることが示唆される．右下にある 3H-7O は静電引力であり，M1+M2 状態を安定化していることが示唆される．
実際に，M1 と M2 の構造上で 3H と 7O の相互作用が確認できるが，N 構造では 3H は 7O ではなく 8O と相互作用していることが分かる．
N+M2 状態と U+M1 状態は左下から右上方向の軸で分離でき，その差異は左下と右上の相互作用を見ればよい．左下の静電引力 1H-10O は N+M2 状態を安定化していることが示唆されるが，実際に N 構造や M2 構造では 1H と 10O が相互作用しているのに対して，U 構造や M1

196

統計数理

第 62 巻

第 2 号 2014

図 9．DIPA 固有値の計算条件依存性．
（a）第 1–10 固有値の刻み幅 Δr 依存性．Δr = 0.05 Å
での固有値に対する相対誤差．rc = 25 Å を用いた．
（b）第 1–10 固有値のカットオフ（a）と（b）では 100 ピコ秒のベリー法を用い距離 rc 依存性．Δr = 0.05 Å を用いた．
た MD シミュレーションの結果を用いた．（c）固有値のベリー法を用いた MD シミュ．Δr = 0.05 Å および rc = 25 Å
レーションの実行時間依存性（Koyama et al., 2011）
を用いた．

構造では離れていることが確認できる．以上のようにして，生体分子内相互作用の摂動解析である PEPCA の結果をバイプロットを用いて表示することで，主成分得点により生体分子の構造状態の同定と，固有ベクトル成分により状態間の差異を生み出している生体分子内相互作用を詳細に明らかにすることができる．また，各状態間に対する安定化あるいは不安定化の寄与
は，静電相互作用が引力的であるかまたは斥力的であるかにより判別できる．
3.2 距離に依存した生体分子・溶媒分子間相互作用の摂動解析（DIPA）の適用
次に，距離に依存した生体分子・溶媒分子間相互作用の摂動解析である DIPA を用いた解析を
（図 6
（a）
）
行う．前節で用いた 1 マイクロ秒のデータから 1 ナノ秒ごとの構造を初期構造としてのみを動かす 500 ピコ秒の MD シミュレーション
ベリー法を用いた溶媒分子（水，Na+ ，Cl− ）（b））．平均溶媒原子数は周期境界の隣接する領域も考慮に入れを合計 n = 1000 本行った（図 6て推定し，fk (r|q) を計算した（図 6（c））．シニョリンと溶媒の置換対称性を考慮した相互作用（d））
は離散化した関数データに対する PCA により行っの数は M = 798 であった．FPCA（図 6．カットオフ距離 rc = 25Å と距離 r の離散化の刻み幅 Δr = 0.05Å を用いた場合にた（付録 B）
は，離散化したデータの変数の数は p = M rc /Δr = 399,000 と大きく，SVD を用いた PCA が実行できなかった．一方で n = 1000 と小さいため，n × n 中心化グラム行列の対角化が可能で（付録 A）．
あり，これを用いて PCA を実行した
まず，DIPA を実装するときに必要となるパラメーターの値を変えたときに固有値がどのよ（a）
は FPCA を離散化した PCA で推定するときに導入した刻み
うに変化するかを調べた．図 9
幅 Δr 依存性である．Δr ≤ 0.2 Å であれば固有値が収束しており，離散化の影響が無くなるこ（b）は式（2.9）で導入したカットオフ距離 rc 依存性である．rc ≥ 20 Å で固有とが分かる．図 9（c）は平均溶媒原子数値が収束しており，カットオフによる影響が無くなることが分かる．図 9
NiJk (r|q) を推定するために必要な溶媒分子のみ動かす MD シミュレーションの実行時間依存性である．100 ピコ秒ほど実行すれば第 1, 2 固有値は収束しており，第 3 固有値もほぼ収束していることが分かる．以上のように，十分な精度のパラメーターを用いればその値に依存しない，DIPA 固有の結果が得られることが分かる．以下の解析では，Δr = 0.05 Å，rc = 25 Å，100 ピ
コ秒の MD シミュレーションでの結果を用いた．

生体分子シミュレーションの摂動解析

197

図 10．DIPA バイプロット（Koyama et al., 2011）．小さい点は主成分得点，曲線は固有関数成分であり，曲線上の数字はシニョリン原子からの距離 [Å] を表す．10O-H などはシ（O）と水の水素原子（H）との静電相互作用をニョリンの 10 番目のアミノ酸の酸素原子表す．第 1, 2, 3 固有関数成分の寄与が上位 5 位までに入る相互作用の一部を示す．ここに示していない相互作用のバイプロットや固有関数成分の寄与を評価する計算法については Koyama et al.（2011）にある．

図9
（b）
（
，c）を見ると，第 1–3 固有値が大きいことから第 1–3 固有関数成分により，シニョリ
ンの構造状態に影響を及ぼす溶媒との相互作用が得られることが期待される．PEPCA の結果（図 10（a）：点）は 4 つの状態を分離できていることが分と同様に，第 1, 2 主成分得点の散布図
かる．また，第 3 主成分得点は N+M1+M2 状態と U 状態を分離していることが分かる（図 10（b）：点）．
次に，各状態の差異を生み出している溶媒との相互作用を明らかにするためにバイプロット
を用いて解析を行う．2.3 節でのバイプロットについての議論は固有ベクトル成分を固有関数成分に置き換えることで DIPA でもそのまま成り立つが，固有関数成分は距離 r をパラメーター（図 10（a））から
として連続的に変化するため，曲線で表される．第 1, 2 主成分のバイプロット
N 状態と M1+M2 状態を識別している相互作用は 8O-O または 8O-H であり，静電引力 8O-H が
M1+M2 状態を安定化していることが示唆される．M1, M2 構造を見ると，8O 原子は溶媒に露出しているのに対して，N 構造ではシニョリン内部に埋もれており，水による 8O 原子の安定化は M1+M2 状態の方が N 状態より大きいことが確認できる．また，固有関数成分の値は 6Å あ
たりで最大となっていることが分かる．距離が離れていくと個々の水との相互作用が弱まっていく効果と，水の個数が増えていく効果があるが，この両者を合わせた効果の構造変化による（b）から静電引力 7O-H が M1+M2 状態に差が，6Å 付近で最大となるためと考えられる．図 10
対して N 状態を安定化していることが示唆されるが，N 構造では 7O 原子が溶媒に露出しているのに対して，M1, M2 構造では 7O 原子はシニョリン内部に埋もれていることが確認できる．

統計数理

198

第 62 巻

第 2 号 2014同様にして，M1 と M2 の相違を生み出している相互作用は 10O-H と 10O-O であり，静電引
（図 10
（a）
）
．実際に M1 構造では 10O は外力 10O-H により M1 が安定されることが示唆されるに露出しているのに対して，M2 構造では内部で相互作用することで，露出が少なくなっている
（図 10
（b）
）を見ると，U 状態と N+M1+M2 状態のことが分かる．第 1,3 主成分のバイプロット
差異を生み出している分子間相互作用は 3OD-H と 3OD-O であり，静電引力 3OD-H により U状態が安定化されていることが分かる．実際に U 状態の構造を見ると 3OD は外に露出しているのに対して，N, M1, M2 の構造では 3OD は露出が少なくなっていることが確認できる．以上のように，DIPA では生体分子の構造変化による特定の原子の周りの溶媒和状態の変化を識別することで構造状態を分離し，その差異を生み出している距離依存的な溶媒との相互作用を明らかにすることが可能である．
4.

おわりに

本稿では生体分子機械を理解するために相互作用の強弱を変える摂動を導入し，その応答性を生体分子構造の分布変化で定量化する摂動解析について紹介を行った．この結果として，生体分子内相互作用に対する摂動解析として PEPCA を導入し，距離に依存した生体分子・溶媒
分子間の相互作用に対する摂動解析として DIPA を導入した．このとき，生体分子構造の分布を大きく変化させるような相互作用の摂動は固有ベクトル/関数，摂動による構造分布変化の大きさは固有値，摂動による構造の比率変化は主成分得点，で得られた．また，これらの結果をバイプロットで表示することにより，構造状態とその差異を生み出している相互作用を明らかにできることを示した．
付A.

録

PCA の異なる実装法

MD シミュレーションで一般的に用いられている原子位置の PCA（Kitao et al., 1991; Garcı́a,や二面角を用いた PCA（dihedral angle PCA; dPCA; Mu et al., 2005;
1992; Amadei et al., 1993）
では，通常は変数の次元 p が数千以下となり，共分散行列の対角化を行うこと
Altis et al., 2007）
が可能であることが多い．これに対して本稿で述べたポテンシャルエネルギー項を用いた PCAや FPCA を推定するための PCA（付録 B）
では数万以上の変数となり共分散行列の対角化を行うことが困難であることが多い．ここでは共分散行列の対角化を用いない PCA の実装法について簡単に述べる．
構造データとして n 個のデータ q1 , . . . , qn が得られたとする．それぞれのデータについてポ
テンシャルエネルギーや離散化した関数データなどの統計量 fi (q) を計算した後，平均を差し引いた（中心化した）
統計量
(A.1)

n
1
f¯i (q) ≡ fi (q) −fi (qj ), i = 1, . . . , pn j=1からなるデータ行列 X̄ij ≡ f¯i (qj ) を考える．共分散行列 C = X̄ X̄ T /n の対角化 C = ri=1 σi2 ui uTiが実行できる場合にはその固有値と固有ベクトルを計算した後，主成分得点は X̄ T ui と計算できる．
X̄ は特異値分解（singular value decomposition; SVD）により特異値 λi > 0，左特異ベクトル
ui ，右特異ベクトル vi を用いて

生体分子シミュレーションの摂動解析

199

表 1．PCA の固有値，固有ベクトル，主成分得点の計算方法．

X̄ =

(A.2)

r


λi ui viT

i=1と表すことができる．このとき共分散行列は
C=

(A.3)

r

λ2i
i=1

n

ui uTiとなり，共分散行列の固有値と固有ベクトルがそれぞれ中心化データ行列の特異値と左特異ベクトルで表されることが分かる．また，主成分得点はX̄ T ui = λi vi(A.4)となり，特異値と右特異ベクトルで表されることが分かる．以上から中心化データ行列の SVD
により PCA の固有値，固有ベクトル，主成分得点が計算でき，バイプロットに必要な全ての情報が得られることが分かる（Jolliﬀe, 2002）．統計解析ソフトの R では prcomp 関数は中心化（R Coreデータ行列の SVD，princomp 関数は共分散行列の対角化により PCA を実行している．
Team, 2013）
また，p が大きく，SVD や共分散行列の対角化ができない場合でも n が小さく中心化グラム
行列 Ḡij = f¯(qi ) · f¯(qj ) の対角化
Ḡ = X̄ T X̄ =

(A.5)

r


λ2i vi viT

i=1が計算できる場合にはその固有値 λ2i と固有ベクトル vi により，共分散行列の固有値 λ2i /n と主成分得点 λi vi が計算でき，共分散行列の固有ベクトルは ui = X̄ vi /λi と計算できることが分か．以上の結果をまとめたのが表 1 である．
る
（Bishop, 2006）
式（A.1）を用いると中心化グラム行列の要素 Ḡij はグラム行列の要素 Gij = f (qi ) · f (qj ) を用
いて
(A.6)

Ḡij = Gij −

n
n
n
1
1 
1 
Glj −
Gim + 2
Glm
n
n m=1
n
l=1

l,m=1

と書き下せることが分かる．このため，データの中心化をあらかじめ行わなくても，グラム行
により中心化グラム行列を計算し，その対角化を実行すれば共分散行列を求めた後，式（A.6）
列の固有値と主成分得点が計算できることが分かる．カーネル PCA（Schölkopf et al., 1998）はカーネル k(q, q  ) = f (q) · f (q  ) を用いることで f (q) を経由することなくグラム行列を計算し，

統計数理

200

第 62 巻

第 2 号 2014

中心化グラム行列の対角化により PCA を実行する．カーネル法の分野では f (q) は特徴空間（feature space）と言われる．カーネル PCA では f (q) を直接計算しないため，ガウシアンカー
ネル k(q, q  ) = exp(−|q − q  |2 /2σ 2 ) のように特徴空間が無限次元の場合も扱えるという利点がある．
B.

FPCA の実装

ここでは FPCA の最も簡単な実装である関数データの離散化を用いた方法を説明する（Jolliﬀe,
2002; Ramsay and Silverman, 2006）．FPCA では関数データ f (r|q) ≡ (f1 (r|q), . . . , fM (r|q))TM  rcに対して線形結合
k=1 0 λk (r)fk (r|q) dr の 分 散 が 最 大 化 に な る よ う な 係 数 λ(r) ≡T(λ1 (r), . . . , λM (r)) を探索する．積分区間 [0, rc ] を NB 個の区間に離散化し，Δr ≡ rc /NB および rl ≡ lΔr とすると
NB
M  rc
M 


(B.1)
λk (r)fk (r|q) dr ≈λk (rl )fk (rl |q)Δr
k=10

k=1 l=1√√のように近似できる．これを M NB 個の係数 λk (rl ) Δr と変数 fk (rl |q) Δr の線形結合とみ√なし，通常の PCA を適用した結果の固有ベクトルを ui (rl ) Δr，固有値を σi2 ，主成分得点を√√TB
gi (q) とする．このとき，PCA の固有ベクトルの直交条件 N
l=1 (u i (rl ) Δr) u j (rl ) Δr = δi,j
は Δr → 0 とすれば，FPCA の固有関数の直交条件 ui , uj  = δi,j に収束することが分かる．また，PCA の主成分得点(B.2)

gi (q) =

NB


√
√
√
(ui (rl ) Δr)T (f (rl |q) Δr − f (rl |q) Δr)l=1

も Δr → 0 とすれば，FPCA の主成分得点 gi (q) = ui (r), f (r|q) − f (r|q)) に収束することが分かる．PCA および FPCA の固有値は主成分得点の分散であることから，Δr → 0 により PCAの固有値も FPCA の固有値に収束することが分かる．以上のように，関数データを離散化後，√f (rl |q) Δr に対して PCA を実行すれば，その固有値と主成分得点がそのまま FPCA の対応す√√る推定値となる．また，PCA の固有ベクトル ui (rl ) Δr を Δr で割ることにより，FPCA の固有関数 ui (r) が推定できる．
謝

辞

本研究は小林徹也博士，上田泰己博士，泰地真弘人博士，友田修司博士，高野光則博士との議論やサポートがなければ遂行することができませんでした．ここに感謝を致します．本稿は日本学術振興会多国間国際研究協力事業「次世代計算機・次世代アルゴリズムによる巨大生体分子システムの動的モデル化」の助成を受けたものです．また，数値シミュレーションの実行は，理化学研究所情報基盤センターの RIKEN Integrated Cluster of Clusters（RICC）システムを用いて行われました．

参

考

文

献

Altis, A., Nguyen, P. H., Hegger, R. and Stock, G.（2007）. Dihedral angle principal component analysis
of molecular dynamics simulations, Journal of Chemical Physics, 126, 244111.
Amadei, A., Linssen, A. B. M. and Berendsen, H. J. C.（1993）. Essential dynamics of proteins, Proteins,17, 412–425.

生体分子シミュレーションの摂動解析

201

Bishop, C. M.（2006）. Pattern Recognition and Machine Learning, Springer, New York.
Case, D. A., Darden, T. A., Cheatham, T. E., III, Simmerling, C. L., Wang, J., Duke, R. E., Luo, R.,
Crowley, M., Walker, R. C., Zhang, W., Merz, K. M., Wang, B., Hayik, S., Roitberg, A., Seabra,G., Kolossváry, I., Wong, K. F., Paesani, F., Vanicek, J., Wu, X., Brozell, S. R., Steinbrecher,T., Gohlke, H., Yang, L., Tan, C., Mongan, J., Hornak, V., Cui, G., Mathews, D. H., Seetin,
M. G., Sagui, C., Babin, V. and Kollman, P. A.（2008）. AMBER 10, University of California,
San Francisco.
Gabriel, K. R.（1971）. The biplot graphic display of matrices with application to principal component
analysis, Biometrika, 58, 453–467.
Garcı́a, A. E.（1992）. Large-amplitude nonlinear motions in proteins, Physical Review Letters, 68,
2696–2699.
Gower, J. C. and Hand, D. J.（1996）. Biplots, Chapman & Hall, London.
Honda, S., Yamasaki, K., Sawada, Y. and Morii, H.（2004）. 10 residue folded peptide designed bysegment statistics, Structure（London）, 12, 1507–1518.
Ikeguchi, M., Ueno, J., Sato, M. and Kidera, A.（2005）. Protein structural change upon ligand binding:Linear response theory, Physical Review Letters, 94, 078102.
Jolliﬀe, I. T.（2002）. Principal Component Analysis, 2nd ed., Springer, New York.
Kitao, A., Hirata, F. and Go, N.（1991）. The eﬀects of solvent on the conformation and the collective
motions of protein: Normal mode analysis and molecular dynamics simulations of melittin inwater and in vacuum, Chemical Physics, 158, 447–472.
古明地勇人, 上林正巳, 長嶋雲兵（2000）. 生体分子の分子動力学シミュレーション（1）
方法, Journal ofChemical Software, 6, 1–36.
Koyama, Y. M., Kobayashi, T. J., Tomoda, S. and Ueda, H. R.（2008）. Perturbational formulation ofprincipal component analysis in molecular dynamics simulation, Physical Review E, 78, 046702.
Koyama, Y. M., Kobayashi, T. J. and Ueda, H. R.（2011）. Perturbation analyses of intermolecularinteractions, Physical Review E, 84, 026704.
Kullback, S.（1959）. Information Theory and Statistics, John Wiley & Sons, New York.
Malolepsza, E., Strodel, B., Khalili, M., Trygubenko, S., Fejer, S. N. and Wales, D. J.（2010）. Symmetrization of the AMBER and CHARMM force ﬁelds, Journal of Computational Chemistry,
31, 1402–1409.
Ming, D. and Wall, M. E.（2005）. Quantifying allosteric eﬀects in proteins, Proteins, 59, 697–707.
Ming, D. and Wall, M. E.（2006）. Interactions in native binding sites cause a large change in protein
dynamics, Journal of Molecular Biology, 358, 213–223.
Mu, Y., Nguyen, P. H. and Stock, G.（2005）. Energy landscape of a small peptide revealed by dihedral
angle principal component analysis, Proteins, 58, 45–52.
R Core Team（2013）. R: A Language and Environment for Statistical Computing, R Foundation for
Statistical Computing, Vienna, Austria, http://www.r-project.org/.
Ramsay, J. O. and Silverman, B. W.（2006）. Functional Data Analysis, 2nd ed., Springer, New York.
Satoh, D., Shimizu, K., Nakamura, S. and Terada, T.（2006）. Folding free-energy landscape of a 10residue mini-protein, chignolin, FEBS Letters, 580, 3422–3426.
Schölkopf, B., Smola, A. and Müller, K.-R.（1998）. Nonlinear component analysis as a kernel eigenvalue
problem, Neural Computation, 10, 1299–1319.
Zwanzig, R.（2001）. Nonequilibrium Statistical Mechanics, Oxford University Press, New York.

202

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

Perturbation Analyses of Biomolecular Simulations
Yohei M. Koyama
Quantitative Biology Center, RIKEN

To understand protein functions in atomic detail, we need to investigate diﬀerent main-chain/side-chain conformational states, interaction patterns to other molecules(compounds, DNA, and other proteins), and contributions of solvents (water and ions).
We review our recent approaches to identify complex conformational states and theirinteractions by performing perturbation analyses of atomic interactions. First, we introduce PEPCA (potential energy principal component analysis), which is a perturbation
analysis of atomic interactions within proteins or complex. PEPCA is implemented by
performing the principal component analysis of interaction energies. PEPCA principal
component scores identify conformational states that have diﬀerent interaction patterns
within proteins or complex. PEPCA eigenvector components identify interactions that
diﬀerentiate each state. Next, we introduce DIPA (distance-dependent intermolecular perturbation analysis), which is a perturbation analysis of atomic interactions between proteins/complex and solvents. DIPA is implemented by performing the functional principal
component analysis by using products of intermolecular forces and average numbers ofsolvent atoms. DIPA principal component scores identify conformational states that havediﬀerent interaction patterns between proteins/complex and solvents. DIPA eigenfunction
components identify interactions that diﬀerentiate each state. Finally, we apply PEPCA
and DIPA to the chignolin folding simulation. Results are visualized by using biplots thatare scatter plots of eigenvector (or eigenfunction) components and principal componentscores.

Key words: Molecular dynamics simulation, perturbation analysis, principal component analysis, functional principal component analysis, biplot.