修士論文 分子アニメーションのための基準振動解析に基づいた 動的

NAIST-IS-MT1351107
修士論文
分子アニメーションのための基準振動解析に基づいた
動的骨格モデルの開発
村岡
伸哉
2015 年 2 月 5 日
奈良先端科学技術大学院大学
情報科学研究科 情報科学専攻
本論文は奈良先端科学技術大学院大学情報科学研究科に
修士(理学)授与の要件として提出した修士論文である。
村岡
伸哉
審査委員：
金谷
重彦 教授
（主指導教員）
佐藤
嘉伸 教授
（副指導教員）
MD.ALTAF-UL-AMIN 准教授
（副指導教員）
杉浦
忠男 准教授
（副指導教員）
佐藤
哲大 助教
（副指導教員）
上野
豊 教授
（産業技術総合研究所）
分子アニメーションのための基準振動解析に基づいた
動的骨格モデルの開発*
村岡
伸哉
内容梗概
分子アニメーションはタンパク質の関与する生化学反応の仮説を、動画とし
て可視化する目的で活用されている。本研究では3Dコンピュータグラフィクス
作成ソフトウェアにおいて、分子モデルを軽処理で動かす骨格の作成手法を開
発する。基準振動を合成した振動運動を骨格に与え、骨格によって分子の熱揺
らぎを表現する。まずタンパク質の原子モデルをその座標と基準振動解析によ
って得られる変位ベクトルを用いて、k-平均法により𝐾個のサブドメインに分割
し、その領域ごとに骨格を決定する手法を考案した。得られた骨格では本来安
定とされる、タンパク質に特徴的な二次構造の1つであるへリックスに歪みが
生じる場合があった。そのため、さらにタンパク質の二次構造を考慮して骨格
を決定する手法に改善した。これらの手法の骨格を設定したヘモグロビンおよ
びアミラーゼのモデルのアニメーションは自然な動きを実現できた。このよう
な骨格モデルを用いた新手法と基準振動解析による原子の変位ベクトルを比較
した結果、骨格モデルによる分子の熱揺らぎの表現は妥当であることが確認さ
れた。
キーワード
3 次元コンピュータグラフィクス、タンパク質、分子アニメーション、分子モデル、基
準振動解析
*奈良先端科学技術大学院大学情報科学研究科情報科学専攻修士論文,
NAIST-ISMT1354407, 2015 年 2 月 5 日.
i
Development of the dynamic skeleton
for molecular mo de ls base d on no rmal mo de
applicable to mole cular mo vie s *
Shinya Muraoka
Abstract
Molecular animation is utilized to visualize the dynamics of chemical reactions
catalyzed by proteins. In this study, we proposed the two methods of generating
skeletons which are embedded in 3D molecular model to move it by thermal
fluctuation. First of all, we developed a method of generating skeletons in every
region in the model. The 𝑘 -means clustering was used to group the atoms into
regions using their coordinates and displacement vectors. The displacement vectors
were calculated from normal mode analysis by using the atomic coordinates of protein.
This method, however, has the problem of distorted motions in -helixes when
performing it till making animations. So next, to overcome this problem, we modify
the method to generate skeletons according to the secondary structure of proteins.
Finally, we confirmed that the models of hemoglobin and amylase behave naturally
when visualizing with the modified method. We also confirmed the displacement
vectors of atoms under both of the proposed methods reasonably correspond with that
of normal mode analysis.
Keywords:
3D computer graphics, protein, molecular animation, molecular model, normal mode
analysis
*Master’s Thesis, Department of Information Systems, Graduate School of Information
Science, Nara Institute of Science and Technology, NAIST-IS-MT1351107, February 5, 2015.
ii
目次
1.
2.
はじめに.......................................................................................................................... 1
1.1
研究の目的............................................................................................................... 2
1.2
提案手法の概要 ....................................................................................................... 2
分子アニメーションの現状 ........................................................................................... 3
2.1
分子アニメーションとその制作環境 ..................................................................... 3
2.2
3DCG ソフトウェア ............................................................................................... 5
2.3
モーションキャプチャ法 ........................................................................................ 6
2.4
分子アニメーション制作を支援するプラグイン .................................................. 7
2.5
分子動力学シミュレーション ................................................................................ 8
2.6
基準振動解析 ........................................................................................................... 9
2.6.1
2.7
3.
k-平均法（k-means 法） ..................................................................................... 11
方法 ............................................................................................................................... 12
3.1
サブドメイン法による骨格の作成 ....................................................................... 12
3.1.1
サブドメインの決定 ...................................................................................... 12
3.1.2
サブドメインの骨格決定............................................................................... 14
3.1.3
振動モードの合成による熱揺らぎの再現 .................................................... 16
3.2
二次構造法による骨格の作成 .............................................................................. 17
3.2.1
DSSP（Define Secondary Structure of Proteins）法と二次構造の分類 17
3.2.2
DSSP による分類に従った二次構造・非二次構造領域の決定 .................. 18
3.2.3
二次構造・非二次構造領域の骨格決定 ........................................................ 19
3.2.4
骨格の平行移動・回転の決定 ....................................................................... 22
3.3
4.
基準振動解析結果のデータベース ............................................................... 10
骨格の定量的評価の方法 ...................................................................................... 23
結果・考察 .................................................................................................................... 25
4.1
サブドメイン法の骨格による結果 ....................................................................... 25
4.1.1
骨格の定性的評価 .......................................................................................... 25
4.1.2
骨格の定量的評価 .......................................................................................... 27
4.2
二次構造法の骨格による結果 .............................................................................. 30
4.2.1
骨格の定性的評価 .......................................................................................... 30
iii
5.
4.3
サブドメイン法と二次構造法の比較 ................................................................... 32
4.4
作成したアニメーションの生物学的な考察 ........................................................ 36
4.5
適用の拡大に関する考察 ...................................................................................... 38
4.6
今後の展望............................................................................................................. 39
おわりに........................................................................................................................ 40
謝辞 ...................................................................................................................................... 41
参考文献 ............................................................................................................................... 42
iv
図目次
図 1. 3 次元座標系と骨格 .............................................................................................. 6
図 2. データベースから取得した基準振動解析結果（最低振動数モード 1） ........ 13
図 3．サブドメインの骨格情報 .................................................................................. 15
図 4．二次構造と非二次構造領域の骨格決定の流れ ................................................ 18
図 5. タンパク質のアミノ酸配列 ............................................................................... 19
図 6. アミノ酸の構造 .................................................................................................. 19
図 7. 3 点位置合わせによる平行移動・回転の計算................................................... 22
図 8．骨格による原子の変位と基準振動による変位の差の計測法 ......................... 23
図 9. サブドメイン法の骨格をヘモグロビンのモデルに設定した結果 ................... 25
図 10. 骨格によるモデルの二次構造変形の仮説 ...................................................... 26
図 11. 振幅の倍率による誤差の変化 .......................................................................... 28
図 12. 最低振動数モードごとの誤差 ......................................................................... 29
図 13．設定された骨格の比較 .................................................................................... 30
図 14. 二次構造法の骨格をアミラーゼのモデルに設定した結果 ............................ 31
図 15. サブドメイン法と二次構造法の骨格による誤差の比較 ............................... 33
図 16. サブドメイン法の骨格（第 2 振動モード）による誤差の色分け ................ 34
図 17. 二次構造法の骨格（第 2 振動モード）による誤差の色分け ........................ 35
表目次
表 1. DSSP による二次構造の分類表 ......................................................................... 17
表 2．ヘモグロビンの最低振動数 1 番目のモードにおける骨格による原子の変位と
基準振動解析結果との変位の誤差 ...................................................................... 27
v
1. はじめに
タンパク質は生物にとって、必須な生体高分子の 1 つである。タンパク質は生体
内において、細胞の構造維持、筋収縮、化合物の輸送、免疫、代謝など、非常に多く
の役割を担っている。タンパク質の機能について知見を広めることが、生物を本質
的に理解するうえで重要となる。そして、タンパク質が働く様子を直接観察するこ
とができれば、機能理解の最たる近道となる。しかし、現在では微小な物体を視覚的
に拡大して観察する顕微鏡を用いても、物理的な限界により、鮮明にタンパク質を
観察することは困難である。そのため、タンパク質の構造から推測される機能の仮
説は計算科学に基づき、時系列ごとの構造変化を予測する分子シミュレーションに
よって、確かめられてきた。一過的に顕れるタンパク質の構造変化から、化合物との
反応機構を考察してきたのである。一方で、DNA の遺伝情報を仲介物質としての
RNA に移す転写や RNA からタンパク質を合成する翻訳、筋収縮など、タンパク質
が関与する連続的な反応の仮説を、分子シミュレーションによって検証することは
困難な課題と考えられ、2 分子のドッキング計算によって個々の分子の相互作用を
評価する計算が活用されている。
こうして、生化学実験や計算機実験によって明らかになった仮説は、3D コンピュ
ータグラフィクスソフトウェア（3DCG ソフトウェア）を用いて作成される分子ア
ニメーションによって、表現されてきた[1]。高速で大容量のメモリを持つ計算機が
低価格化し一般に利用できるとは言え、アニメーションを作成する際に、膨大な計
算量と計算時間を要する分子シミュレーションを 3DCG ソフトウェア上で直接適用
することや、大量の計算結果をそのまま利用することは難しかった。
そのような背景のなか、分子アニメーションでは祖視化したモデルを導入して、
模式的にモデルを動かす表現が有効な場合もあり、軽処理で科学的な予測に基づき
モデルを動かすことが要求されてきた。
1
1.1 研究の目的
これまで、3DCG ソフトウェアにおいては、人体のモデルを動かす骨格が利用さ
れてきた。本研究では、これまで分子モデルには利用されていなかった骨格による
動画作成手法を開発した。この骨格には、基準振動解析の結果から得られる分子の
振動モードを合成した熱揺らぎの運動を動きとして与える。分子モデルを計算量の
多いシミュレーションではなく、骨格によって生体内で予想される動きで動かすこ
とを目的とする。
1.2 提案手法の概要
本提案手法の概念は“原子の動きを粗視化するために平均化した骨格を作成し、
それによって軽処理で分子モデルを動かす”ということである。本研究では 3DCG
ソフトで利用可能な骨格を作成し、基準振動解析の結果を用いて動きを与える ため
の 2 つの手法を提案する。この骨格オブジェクトを 3DCG ソフトウェア上に表示さ
せた分子モデルに設定することで、骨格の動きに追従させて分子モデルを動かすこ
とが可能となる。これにより、分子モデルを計算科学に基づいて動かすことができ、
3DCG ソフトウェア上における計算量を減らすことで、多数の分子モデルを扱う分
子アニメーション制作の円滑化を実現できる。
基準振動解析の結果には、分子内の全原子のそれぞれの変位ベクトルが含まれる。
説明する 2 つの提案手法では、いくつかの原子の変位ベクトルを平均化して、それ
を骨格の動きとし分子モデルを動かすものである。
2
2.
分子アニメーションの現状
まず現在の分子アニメーションとその制作に用いられるソフトウェアについて説
明する。次にソフトウェアの機能を支える技術や理論について、説明を加える。
最後に現在制作で生じている課題に対処するため、本研究で提案する 3DCG ソフ
トウェアの機能を使った、タンパク質分子のモデルを動かす骨格を作成する新規手
法の概要について説明する。また、その際に必要となる技術についても述べる。
2.1 分子アニメーションとその制作環境
生体内高分子の関与する化学反応を高機能な 3DCG ソフトウェアを用いて、動画
として可視化する分子アニメーションは、現実には観測できないタンパク質の動き
を直感的に表現することができ、時間的に変化する酵素の立体構造とそれに触媒さ
れる化学反応を分かりやすく表現できるという利点がある。そのため分子アニメー
ションを用いると、初学者や研究者が分子世界を容易に理解することができる。ま
た、生物学を専門に学んだことの無い一般の人々に、研究内容を伝えるためのツー
ルの一つにも成り得る。しかし、研究者にとって、煩雑な操作が要求される 3DCG
ソフトウェアを駆使して分子アニメーション制作を行うのは非常に手間のかかる作
業である。
それでも、必要な機能を備えたハードウェアの低価格化により、計算機を使った
分子アニメーション制作に対する障壁は低下した。また、ソフトウェアの拡張機能
を用いて、制作を円滑化するための 3DCG ソフトウェア用の機能を追加するプラグ
インや、粗視化された分子モデル表示と制作に必要最低限な機能が整備されたソフ
トウェア Molecular Flipbook 1も開発されている。また数種の化学反応のアニメーシ
ョン制作について、チュートリアルがまとめられたサイト Molecular Movies2も存
在し、制作は活発になりつつある。
しかし、分子アニメーション制作に用いられるタンパク質分子のモデルは、タン
パク質構造データバンク 3が規定する PDB ファイルと呼ばれるファイルに格納され
た静的な 3 次元の構造情報フォーマットデータ（Protein Data Bank file format data）
1 https://www.molecularflipbook.org
2 http://www.molecularmovies.com
3 http://www.rcsb.org/pdb/home/home.do
3
を基に作製されるため、その動きは実証されていない仮説に基づいて加えるのが現
状である。また、3DCG ソフトウェアのプラグインには分子の動きを MD シミュレ
ーションによって求める機能もあるが、正常に計算が行われない等の不具合が解消
されず、利用者による対処も難しいのが現状である。さらに、アニメーション制作で
用いられるタンパク質の分子数を増加した系を扱うと、MD 計算を適用させる計算
量が膨大となるために一般的な制作環境では難しい。
このような背景のなかで、生体内で予測される動きに近づけつつも、軽処理で分
子モデルを動かす骨格を用いた方法はまだ試みられていなかった。
4
2.2 3DCG ソフトウェア
3DCG ソフトウェアには Maya 4, Cinema4D 5, Shade6等の商用ソフトだけでなく、
フリーソフトである Blender 7も高機能で広く利用されている。3DCG ソフトウェア
における立体構造の表現は幾何学的モデルであり、ポリゴンメッシュによる形状デ
ータとテクスチャ画像の貼付によって構築される。また、柔らかい物体や液体の自
然な動作を表現できる物理シミュレーションを元にしたパーティクル法の利用も普
及している[2]。
Blender においては機能拡張のためのプラグインを個人が作成することも可能で
あり、作成された拡張機能は、公開プログラム開発支援サイト GitHub 8において公
開されている。これらはオープンソースとして改変や機能追加が可能であり、開発
者コミュニティによっては、さかんに改善や開発が進められ、商用ソフトウェアに
比べて開発スピードが早いのが特徴である。また、分子アニメーション制作を支援
する PDB データを利用できる拡張モジュール（ePMV9[3], BioBlender10[4]）も存在
する。
一般に 3DCG ソフトウェアはポリゴンによる立体モデルを基本とするため、原子
モデルの扱いは必ずしも充分ではなく、原子数が多いモデルでは難しい。そのため、
タンパク質を構成する一部の原子を模式的に表示するなどして、対処する必要性が
生じる。また、プラグインを用いず 3DCG ソフトウェアで分子モデルを利用する場
合は、タンパク質や低分子の立体構造情報が保存された PDB データを分子モデリン
グソフトウェアで表示させ、3DCG ソフトウェアで扱えるフォーマットに変換する
必要がある。
4 http://www.autodesk.co.jp/products/maya/overview
5 http://www.aanda.co.jp/products/cinema4d/
6 http://shade3d.jp/
7 http://blender.jp/
8 http://www.bioblender.eu
9 http://epmv.scripps.edu
10 https://github.com
5
2.3 モーションキャプチャ法
人物に貼り付けた骨格マーカーの動きをビデオ撮影し、その動きを骨格データと
して記録する手法はモーションキャプチャと呼ばれる。この記録データには、人体
モデルに適用できる骨格のアニメーションにおけるフレームごとの位置と角度が記
述されており、3DCG ソフトウェアで利用できるようになっている。人間のリアル
な動きが保存された骨格は、人体のモデルに設定され、映画などのコンピューター
ないメーションやゲームに登場するキャラクターの人間らしい動きの再現に利用さ
れる。3DCG ソフトウェアにおいて、モデルのポリゴンメッシュへの骨格の関連付
けによって、リアルタイム変形を行う処理はスキニングと呼ばれる 11。ポリゴンの頂
点を骨格に割り当てることにより、骨格の動きに追従させポリゴンを移動する。複
数の骨格の寄与を重み付けし、スムーズな動きを実現する手法が利用されている[5]。
BVH ファイルはモーションキャプチャデータを格納するための代表的なファイ
ルフォーマットであり、BioVision 社 12によって策定された。BVH ファイルは図 1 に
示すような 3 次元座標系における 1 つの骨格の、1 フレーム単位での動作を記述し
たものである[6]。フレームごとに骨格の位置と回転を記述することで、3DCG ソフ
ト側でフレーム間の値は補完（3DCG ソフトによっては補完の設定が変更可能 13）さ
れるので、骨格は滑らかな動きを表現できる。
図 1. 3 次元座標系と骨格
11 http://www.oakcorp.net/polytrans/boneskel.php
12 http://www.biovision.com
13 http://wiki.blender.org/index.php/Doc:JA/2.6/Manual/Animation/Editors/Graph/ FCurves
6
2.4 分子アニメーション制作を支援するプラグイン
ePMV14は、3DCG ソフトである Cinema4D, Maya や Blender のためのフリーの
プラグ イン である 。3DCG ソフト ウェ アで 分 子モデ リン グソフ トウ ェア であ る
Python Molecular Viewer[7]を直接実行することが可能であり、これにより 3DCG
ソフトウェアのインターフェース上で、PDB データを用いた分子モデルを表示させ
ることができる。ePMV が他の分子ビューワーと比較して優れている点は、ホスト
である 3DCG ソフトウェアの機能によって、読み込んだモデルの色や物体表面の質
感の変更、シェーディングなどが可能な点である。さらにそれらの情報をレンダリ
ングすることによって高精細な画像や動画を作成できる。しかし、PDB データを読
み込んでモデルを表示する際や、モデル自体の表示方法（リボンモデル、分子表面モ
デル、ボール＆スティックモデルなど）を変更する際に相応の計算時間を要する。ま
た、ePMV には分子モデルを計算によって動かす機能もあるが、分子ビューワーと
比較して十分に機能しないという問題もある。また、表示させるモデルの原子数が
増加することにより 3DCG ソフト自体の処理が低下し、モデリング機能も鈍化する。
14 ePMV, http://epmv.scripps.edu
7
2.5 分子動力学シミュレーション
分子動力学シミュレーションにより、タンパク質やDNA などの生体高分子の構
造・物性・機能相関を情報学的に解明することができる[8]。分子動力学シミュレー
ションに良く用いられるアルゴリズムは古典分子動力学法（古典Molecular
Dynamics、古典MD）である。古典MDは、 一般的に「N 体多体系の座標の時間
経過を、ニュートンの運動方程式の数値解を差分法により求めてシミュレーション
すること」と定義される。
古典MDの大まかな手順は以下に示す。古典MDの計算には、分子の初期構造と
力場定数を用いる。初期構造（座標）は、Protein Data Bank（PDB） 15などに登
録されているX 線結晶解析や多次元NMR による3次元座標が用いることが多い。
ただし、PDBの座標を使用する場合は、足りない原子（特に水素）や残基を足した
り、イオンを加えたりする必要がある。
(1). 初期構造（座標）を入力し、速度を発生させる。
(2). その座標と力場定数を用いて力を計算する。
(3). 時間刻み（time step, 𝛿𝑡 ） 後の座標と速度を求める。
(4). 必要に応じて温度や圧力の制御を行う。
（2）～（4） の各段階を繰り返すことにより、各原子の座標と速度、さらに座
標と速度から得られるエネルギーなどの物理量が時系列で得られる。これらを用い
て処理や画像解析等を行うことにより、対象とする生体高分子の構造や物性につい
て議論することができる。なお、各原子の三次元座標の時間変化はトラジェクトリ
（trajectory, 軌跡）と呼ばれる。MD計算は非常に計算量が多く、数千原子で構成
されるタンパク質単体の定常状態における動作を数p（ピコ、10-12）秒解析するの
にも数日以上の計算時間がかるのが通常である。さらに複合体の相互作用において
は、相互作用する相互配置等の計算条件を仮定する必要があり、動きを観察しても
得られるものは殆どなく、仮説の有効性を検証するのに利用されるのが一般的であ
る。
15 http://www.rcsb.org/pdb/home/home.do
8
2.6 基準振動解析
生体分子のダイナミクス解析のための手法の一つである基準振動解析は生体分子
もミクロな「構造物」とみなし、その振動モード（振動の形）を工学的に解析するこ
とで動特性を知ろうとするものである。生体分子の運動は複雑であるが、基準振動
解析を用いることで、その運動の特徴を少数のモードで記述することができる[9]。
弾性ネットワークモデルによる基準振動解析では、各原子の運動を近似的に調和
振動子（基準振動）の和で記述し、エネルギー極小点近傍での小さな振動を対象と
する。分子内振動の各自由度𝑛に対応した座標（𝑞1 , 𝑞2 , … , 𝑞𝑛 ）、ベクトル表示を𝒒と
し、エネルギー極小点近傍でエネルギー関数E(𝒒)が多次元放物面で近似できると考
えると、原子𝑖と𝑗の間のポテンシャル関数は次のような単純な形で与えられる
[10]。
−(𝑑𝑖𝑗 )2
𝑃𝐷𝐵 2
} (𝑑𝑖𝑗 − 𝑑𝑖𝑗
)
E(𝑑𝑖𝑗 ) = 𝐾𝑒𝑥𝑝 {
2
𝑎
(1)
𝑃𝐷𝐵
𝑑𝑖𝑗 は計算している分子構造における原子𝑖と𝑗の間の距離である。𝑑𝑖𝑗
はPDBにおけ
る原子𝑖と𝑗の間の距離である。𝐾と𝑎は原子の種類によらないパラメータである。安
定点での二次微分行列（Hessian）𝐹𝑖𝑗 はエネルギーから算出され、運動方程式に利
用される。
𝜕2𝐸
𝐹𝑖𝑗 =
𝜕𝑞𝑖 𝜕𝑞𝑗
(2)
基準振動解析では、原子座標のゆらぎが調和振動子であって、このエネルギー表面
においても存在すると仮定すると、運動方程式は定係数の線形微分方程式であるた
め、その一般解𝑞𝑖 は固有振動の重ね合わせとして、係数𝐴𝑖𝜆 と固有振動数𝜔𝜆 を用いて
次式のように表すことができる。
𝑁
𝑞𝑖 =
𝑞𝑖0
+ ∑ 𝐴𝑖𝜆 𝑄𝜆
𝜆=𝑙
𝑄𝜆 = 𝛼𝜆 sin(𝜔𝜆 𝑡 + 𝛽𝜆 )
9
(3)
式（3）中の𝑄𝜆 は基準振動である。𝛼𝜆 は振幅、𝜔𝜆 は固有振動数、𝛽𝜆 は位相のずれを
表す。
式（3）は、座標𝑞𝑖 には個々の原子運動が直接対応しているのに対して、基準振
動𝑄𝜆 はタンパク質分子内の全原子が示す集団運動に対応する。行列𝐹の固有値を求
めることで、𝜔𝑖2 が固有値として得られ、対応する固有ベクトルが重ね合わせの係数
𝐴𝑖𝜆 となる。得られた𝜔𝜆 と𝐴𝑖𝜆 を式（3）に代入することで生体内高分子の熱揺らぎ
を示す原子の運動が計算できる。式（1）によるポテンシャル関数は、分子力場に
基礎を置いて定義されたものではなく、統計熱力学的に正確な温度をとらないこと
から、𝛼を決めることはできないので、系の全原子の振幅の平均を一定にする値を
用いる[10]。
2.6.1
基準振動解析結果のデータベース
立体構造解析によって明らかになった約 10 万種のタンパク質の構造に対しての
基準振動解析の結果が Promode Elastic 16[11]と呼ばれるデータベースに蓄積されて
いる。基準振動解析の計算には、輪湖らの開発したプログラムである PDBETA
（Elastic-network-model based normal mode analysis in Torsional Angle space
for PDB data）が使用されており、このプログラム自体もデータベースから取得可
能である。データベースに公開されていないタンパク質についても個人で解析を行
うことが可能であるが、原子数 4000 を越えるタンパク質になると計算量が増大す
る。計算時間は数分程度であり、一般的な PC 上でも容易に動かすことができる。
Promode Elastic データベースのそれぞれのタンパク質のページでは、分子ビュ
ーアーである jV17や Jmol18を用いたアニメーションや原子の変位ベクトル図を確認
できる。1 分子単体の基準振動解析結果を含めた動的な振る舞いが観察可能である。
16 http://promode.pdbj.org/promode_elastic/index.do
17 http://pdbj.org/jv
18 http://jmol.sourceforge.net
10
2.7 k-平均法（k-means 法）
k-平均法の目的は d 次元の𝑁個のデータ𝐷 = {𝑥1 , … , 𝑥𝑁 }, 𝑥𝑖 ∈ 𝑅𝑑 をデータ間の類似性
（距離）を尺度に、あらかじめ定めた𝐾個のクラスタに分類することである[12]。k平均法のアルゴリズムを以下に示す。
初期化： 𝑁個のデータをランダムに𝐾個のクラスタに振り分け、それぞれのクラス
タの平均ベクトル（重心、セントロイド）を求め、𝝁𝑘 (𝑘 = 1, … , 𝐾)とする。
（1）i 番目のデータが、k 番目の代表ベクトルが支配するクラスタに帰属するか否
かを表す帰属変数𝑞𝑖𝑘 を次式のように決定する。
𝑞𝑖𝑘
||𝑥𝑖 − 𝜇𝑗 ||2 の場合)
1 (𝑘 = arg min
𝑗
={
0
(それ以外の場合)
(4)
（2）次式にしたがって、セントロイド𝝁𝑘 を決める。
𝑁
∑𝑁
𝜕𝐽(𝑞𝑖𝑘 , 𝝁𝑘 )
𝑖=1 𝑞𝑖𝑘 𝒙𝑖
= 2 ∑ 𝑞𝑖𝑘 (𝑥𝑖 − 𝝁𝑘 ) = 0 ⇒ 𝝁𝑘 = 𝑁
∑𝑖=1 𝑞𝑖𝑘
𝜕𝑥
(5)
𝑖=1
（3）繰り返し：上記（1）と（2）を次式の評価関数が収束するまで（状態変化がな
くなるまで）繰り返す。
𝑁
𝐾
𝐽(𝑞𝑖𝑘 , 𝝁𝑘 ) = ∑ ∑ 𝑞𝑖𝑘 ||𝑥𝑖 − 𝝁𝑘 ||2
(6)
𝑖=1 𝑘=1
この最適化アルゴリズムの収束先は初期値に依存するので、最適解に近い解を得
るためには何回か初期値を変えて実行する必要がある。
11
3. 方法
この章では 3DCG ソフトウェアで利用されるタンパク質分子モデルに設定する骨
格（BVH ファイル）を作成する 2 つの手法について説明する。
3.1 節でサブドメイン法について、3.2 節で二次構造法について説明を行う。
3.1 サブドメイン法による骨格の作成
タンパク質は数千から数万の原子で構成される。また、基準振動解析によって、分
子がある振動数で振動するときの、それぞれの原子の変位ベクトルが得られる。タ
ンパク質を構成する原子を座標と変位ベクトルを用いて、プログラミング言語 R
（GNU-R） 19のパッケージに含まれる k-平均法によって、 𝐾個のサブドメイン（ク
ラスタ）に分割する。そして、分割したサブドメインごとに骨格を決定する。本研究
ではこの領域をサブドメインと呼ぶので、この骨格を作成する手法をサブドメイン
法と命名した。タンパク質にはドメインと呼ばれる、特定の機能を持った領域も存
在するが、さらにそれを細分化する手法である。骨格の動きは、基準振動解析の結果
を元に最低振動モードの合成によって、平行移動のみで与えられる。
3.1.1
サブドメインの決定
ProMode Elastic データベースに保存されている解析結果のうち、10 個の最低振
動数モードについて、それぞれの原子の振動が記されているファイル
（xxxx_1.flcatmE, xxxx_2.flcatmE, ..., xxxx_10.flcatmE）を使用した。これは同デ
ータベースから取得可能な PDBETA（v1.01）プログラムによっても、パラメータな
どを変化させて計算可能である。本研究ではヘモグロビン（PDB ID: 1A00）の結果
を用いた。結果にはタンパク質の原子の 3 次元空間における座標と変位ベクトルが
含まれており（図 2; Atomic coordinate, Displacement vector）、この 2 つのデータ
を用いて k−平均法を適用させることで全原子を𝐾個のサブドメインに分割した。本
研究では𝐾=16 とした。
19 http://www.r-project.org/
12
図 2. データベースから取得した基準振動解析結果（最低振動数モード 1）
Atomic coordinate（原子座標）と Displacement vector（変位ベクトル）を用い、
近隣で同時に動く原子が同じクラスタ（サブドメイン）に含まれることを期待し、
k-平均法を実行した。
13
使用するデータのスケーリングは（7）式、または（8）式で行った。i 番目の原子
の座標を 𝒓𝑖 (𝑥𝑖 , 𝑦𝑖 , 𝑧𝑖 )とし、m 番目の最低振動モードにおける i 番目の原子の変位ベク
トルを 𝒅𝑚𝑖 (𝑥𝑚𝑖 , 𝑦𝑚𝑖 , 𝑧𝑚𝑖 )とする。このとき、スケーリング後の i 番目の原子の座標 𝒓𝑖 は
𝒁𝑟𝑖 とし（7）式で、m 番目のモードの i 番目の原子の変位ベクトル 𝒅𝑚𝑖 を𝒁𝑑𝑚𝑖 とし（8）
式に示す。（7）式における𝜎𝑥2 , 𝜎𝑦2 , 𝜎𝑧2はタンパク質に含まれる全原子の、x,y,z 成分そ
れぞれの分散である。
（8）式における𝜎𝑥2𝑚 , 𝜎𝑦2𝑚 , 𝜎𝑧2𝑚 は m 番目のモードの全原子の変位
ベクトルの、x,y,z 成分それぞれの分散を表す。
𝒁𝑟𝑖 =
Sr =
𝒁𝑑𝑚𝑖 =
Sr
=(
√𝜎𝑥2
𝒅𝑚𝑖
S𝑑𝑚
S𝑑𝑚 =
3.1.2
𝒓𝑖
𝑥𝑖 𝑦𝑖 𝑧𝑖
, , )
Sr Sr Sr
+ 𝜎𝑦2 + 𝜎𝑧2
(7)
3
=(
𝑥𝑚𝑖 𝑦𝑚𝑖 𝑧𝑚𝑖
,
,
)
S𝑑𝑚 S𝑑𝑚 S𝑑𝑚
√𝜎𝑥2𝑚 + 𝜎𝑦2𝑚 + 𝜎𝑧2𝑚
(8)
3
サブドメインの骨格決定
サブドメインごとの骨格の位置は基準振動解析の結果に含まれる原子座標、変位
は変位ベクトルをもとに決定した。各骨格の変位は最低振動数のモードの個数あ
り、本研究では 10 個のモードの結果を用いているので、骨格ごとに 10 個の変位が
与えられる。骨格の情報については、図 3 に示した。また、骨格の中心と主軸の情
報を用いて初期位置を決定し、BVH ファイルのフォーマットにしたがって記述し
た。すべての変位について合成し、アニメーションにおけるフレームごとに時系列
の変位を求め、同ファイルに記述して骨格の動きとした。変位の合成方法について
は、3.1.3 節にて説明を行う。
14
Z
(b)
(a)
(b)
Y
(d)
X
骨格における名称
a
中心
b
主軸
c
計算式
𝑘
∑𝑁
(
)
𝑖=1 𝒓𝑖 𝑥𝑖 , 𝑦𝑖 , 𝑧𝑖
𝑁𝑘
(𝜎𝑥𝑘 ,
𝑣𝑥𝑘 𝑦𝑘 𝑣𝑥𝑘 𝑦𝑘
,
)
𝜎𝑥𝑘 𝜎𝑥𝑘
振動モード m における
𝑘
∑𝑁
𝑚=1 𝑑𝑚𝑘 (𝑥𝑚𝑘 , 𝑦𝑚𝑘 , 𝑧𝑚𝑘 )
中心からの変位
𝑁𝑘
図 3．サブドメインの骨格情報
（a）骨格の中心。あるクラスタ k に含まれる原子座標の中心（平均）を表す。こ
こで、𝑁𝑘 はクラスタ k に含まれる原子の総数を表す。（b）骨格の始点。中心から
（b）骨格の大きさと向きを決定する軸。クラスタ k における原子の x 座標の標準
偏差𝜎𝑥𝑘 を軸の x 座標とし、xy の共分散𝑣𝑥𝑘𝑦𝑘 、xz の共分散𝑣𝑥𝑘𝑦𝑘 をそれぞれ原子の
x 座標の標準偏差𝜎𝑥𝑘 で除した値を軸の y 座標、z 座標とすることで軸を表す。
（c）
骨格の中心からの変位。クラスタ k における、m 番目のモードの原子の変位ベク
トルの平均。
15
3.1.3
振動モードの合成による熱揺らぎの再現
サブドメインごとの骨格の変位は振動モードの数だけ計算され、全てのモードの
変位を合成することによって熱揺らぎを再現する。k 番目の骨格の時刻 t における
変位は、𝑀番目までのモードの変位を合成することによって、
（9）式によって求めら
れる。今回は最低振動数のモード 10 個までの解析結果から骨格の変位を求めている
ので、𝑀 = 10である。
𝑀
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
𝑊𝑘 (𝑡) = ∑ ⃗⃗⃗⃗⃗⃗⃗
𝑉𝑚𝑘 ∙ 𝐴𝑚 ∙ sin(2𝜋𝜔𝑚 𝑡 + 𝛼𝑚 )
(9)
𝑚=1
⃗⃗⃗⃗⃗⃗⃗
ここで𝑉
𝑚𝑘 は k 番目の骨格の m 番目のモードの変位、𝐴𝑚 は m 番目のモードの振
幅、𝜔𝑚 は m 番目のモードの振動数、𝛼𝑚 は m 番目のモードの位相である。本研究で
は位相は 0、振幅は 5 とした。
16
3.2 二次構造法による骨格の作成
タンパク質の構造情報データである PDB データに記述されている原子座標より、
DSSP 20と呼ばれるアルゴリズムによって、分子内における二次構造が決定される。
この二次構造はへリックスやシートなどの、タンパク質の部分的な立体構造のこ
とである。これらの二次構造を保持するように、二次構造それぞれに 1 つ骨格を割
り当てる手法が二次構造法である。ただし二次構造をつくらない原子は、サブドメ
イン法と同様に原子をその座標と基準振動解析の変位ベクトルにより分類し、サブ
ドメインごとに骨格を決定する。骨格の動きは、平行移動だけでなく回転を含んだ
動きを計算する。
3.2.1
DSSP（Define Secondary Structure of Proteins）法と二次構造の分類
タンパク質の二次構造は原子構造に基づいて決められ、プログラム DSSP として
知られている[13]。基準振動解析データベースである ProMode Elastic の解析結果
においても、その分類が記されている。分類の種類については、表 1 に示した。ま
た、二次構造の詳細な割り当ての方法については、付録に記載した。
表 1. DSSP による二次構造の分類表
Code
Explanation
H
-helix
B
residue in isolated -bridge
E
extended strand, participates in -ladder
G
3-helix (3 10 helix)
I
5-helix (-helix)
T
hydrogen bonded turn
S
bend
(Blank)
strands for loop or irregular
20 http://swift.cmbi.ru.nl/gv/dssp/
17
3.2.2
DSSP による分類に従った二次構造・非二次構造領域の決定
DSSP 21による分類にしたがって、二次構造・非二次構造領域の決定を行う。DSSP
の分類表は、 表 1 に記載した。本研究では H,B,E,G,I の領域を二次構造として取り
扱い、その他の領域については非二次構造領域として扱った。この手法では、骨格に
は平行移動だけでなく回転の動きも与えるため、その計算に必要な座標と変位ベク
トルの算出と抽出も行った。これについては、3.2.3 節で説明する。作業フローを 図
4 に示す。
図 4．二次構造と非二次構造領域の骨格決定の流れ
21 http://swift.cmbi.ru.nl/gv/dssp/
18
3.2.3
二次構造・非二次構造領域の骨格決定
まず、今後の手順の説明に必要な炭素、N 末端、C 末端について説明を行う。タ
ンパク質は図 5 に示すように複数のアミノ酸が結合し、折り畳まることで構造をな
す。ここで炭素とは、タンパク質の主鎖（-NH-CH-CO-）において、赤枠で示すペ
プチド結合（-CO-NH-）に関与しない、青矢印で示す炭素（C）のことを言う。図 6
には 1 つのアミノ酸の構造を示した。また N 末端はアミノ基で、C 末端はカルボキ
シル基で終結している側の末端である。アミノ酸配列を書くとき N 末端は左に置き、
C 末端にかけて書く（図 5）。図 5、図 6 ともにアミノ酸の側差は R で表した。
ペプチド結合
N 末端
C 末端
炭素
炭素
図 5. タンパク質のアミノ酸配列
カルボキシル基
アミノ基
炭素
側鎖
図 6. アミノ酸の構造
19
炭素
二次構造 1 つには 1 つの骨格、非二次構造領域には k-平均法によって分割したサ
ブドメイン 1 つに対して骨格が 1 つ決定されるようにする。
二次構造はそのアミノ酸の個数によって、平行移動・回転の計算に必要な 3 点の
座標と変位ベクトルの計算方法が異なるため、その条件と手順を以下に記した。ア
ミノ酸数（炭素）が 3 個に満たない二次構造については骨格を定義しない。
（1）アミノ酸数が 3 個以上の二次構造を抽出。
（2）以下の分類にしたがい、骨格の平行移動・回転の計算に必要な座標と変位ベク
トルそれぞれの平均の計算、もしくは抽出を行う。
・アミノ酸数≧7
第 1 点：N 末端から 3 アミノ酸の炭素の平均
第 2 点：二次構造の中心となるアミノ酸の炭素
第 3 点：C 末から 3 アミノ酸の炭素の平均
・アミノ酸数 6
第 1 点：N 末端から 1,2 番目のアミノ酸の炭素の平均
第 2 点：N 末端から 3 番目の炭素
第 3 点：N 末端から 5,6 番目のアミノ酸の炭素の平均
・アミノ酸数 5
第 1 点：N 末端から 1,2 番目のアミノ酸の炭素の平均
第 2 点：N 末端から 3 番目の炭素
第 3 点：N 末端から 4,5 番目のアミノ酸の炭素の平均
・アミノ酸数 4
第 1 点：N 末端から 1,2 番目のアミノ酸
第 2 点：N 末端から 3 番目の炭素
第 3 点：N 末端から 4 番目の炭素の平均
・アミノ酸数 3
第 1 点：N 末端から 1 番目の炭素
第 2 点：N 末端から 2 番目の炭素
第 3 点：N 末端から 3 番目の炭素
20
また非二次構造領域の 3 点の座標と変位ベクトルの抽出方法については、以下の
手順にしたがった。クラスタリングによって、1 つのアミノ酸の原子が複数のサブド
メインに含まれていることが考えられたので、サブドメインにアミノ酸の炭素が含
まれているかで、あるアミノ酸がどのサブドメインに属すかを決定した。
（1）3.1.1 節と同様の方法で、座標と変位ベクトルのスケーリングを行う。スケー
リングされたデータを用いて、k-平均法により非二次構造領域の原子を𝐾個のサブド
メインに分割する（𝐾＝40）。
（2）3 点の計算はサブドメインに含まれるアミノ酸が 3 つ以上連続していた場合、
それらを抽出する。また、サブドメインにアミノ酸が 3 つ以上連続する領域が 2 つ
以上あった場合、それらを分割し、抽出する。例えば、サブドメインに N 末からの
アミノ酸番号で 1,2,3 と 6,7,8 のアミノ酸が属していた場合、1,2,3 と 6,7,8 として分
割する。
（3）骨格の平行移動・回転の計算に必要な 3 点の座標と変位ベクトルの抽出、計算
は、前頁の二次構造の分類と同様に、
（2）において抽出したアミノ酸の炭素をもと
に行った。
以上より得られた 3 点の情報は 3.2.4 節において、骨格の平行移動・回転計算に利
用される。
21
3.2.4
骨格の平行移動・回転の決定
3.1.3 節と同様の方法により、骨格の時系列変位が求まる。この変位ごとに、骨格
の平行移動と回転を計算し、骨格の動きとする。回転と平行移動の計算は、3 点位置
合わせを用いる（対応する点は既知）。本研究では基準振動によって原子が大きく変
位することなく、また回転移動量が微小であるから、より簡単な 3 点の位置合わせ
を用いた。その概略図を 図 7 に示す。
図 7. 3 点位置合わせによる平行移動・回転の計算
i,ii,iii をそれぞれ、1,2,3 点目の座標とする。（a）変位前の i が、変位後の i に重
なるように平行移動させる。このときの、距離を求める。
（b）ii と i がすべて一
直線上に並ぶように回転させ、そのときの回転を求める。
（c）iii が最も近くなる
ように、i,ii を並べた線の周りに回転し、そのときの回転を求める。（a）から
（c）の計算によって、骨格の平行移動と回転が求まる。
22
3.3 骨格の定量的評価の方法
骨格による振動が基準振動解析の結果を十分に表現しているかを評価するために、
骨格によって動かされる原子の変位と基準振動解析の結果の変位との比較を行った。
その方法を 図 8 に示す。ここでの m 番目の振動モードにおける、i 番目の原子の誤
差𝑒𝑚𝑖 (𝑡)については（10）式を用いて計算した。m 番目の振動モードで動く骨格によ
る i 番目の原子の時刻 t における変位座標は 𝒓𝑏𝑚𝑖 (𝑡) = (𝑥𝑏𝑚𝑖 (𝑡), 𝑦𝑏𝑚𝑖 (𝑡), 𝑧𝑏𝑚𝑖 (𝑡))で、m 番
目 の 振 動 モ ー ド で 動 く i 番 目 の 原 子 の 時 刻 t に お け る 変 位 座 標 を 𝒅𝒎𝒊 (𝑡) =
(𝑥𝑚𝑖 (𝑡), 𝑦𝑚𝑖 (𝑡), 𝑧𝑚𝑖 (𝑡))で表す。
骨格による原子の変位
𝒓𝑏𝑚𝑖 (𝑡) = 𝑥𝑏𝑚𝑖 (𝑡), 𝑦𝑏𝑚𝑖 (𝑡), 𝑧𝑏𝑚𝑖 (𝑡)
𝒆𝒎𝒏 (𝑡)
基準座標
基準振動解析の結果の変位
（PDB データの座標）
𝒅𝒎𝒊 (𝑡) = (𝑥𝑚𝑖 (𝑡), 𝑦𝑚𝑖 (𝑡), 𝑧𝑚𝑖 (𝑡))
図 8．骨格による原子の変位と基準振動による変位の差の計測法
𝑏
𝑏
𝑏
𝑒𝑚𝑖 (𝑡) = √(𝑥𝑚𝑖
(𝑡) − 𝑥𝑚𝑖 (𝑡))2 + (𝑦𝑚𝑖
(𝑡) − 𝑦𝑚𝑖 (𝑡))2 + (𝑧𝑚𝑖
(𝑡) − 𝑧𝑚𝑖 (𝑡))2
23
(10)
今回は作成した骨格に、最低振動数第 1 番目のモードの動きを与え、検証を行っ
た。3.1.3 節の合成式より、k 番目の骨格の 1 番目のモードの動きは次式によって得
られる。
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
𝑊𝑘 (𝑡) = ⃗⃗⃗⃗⃗⃗
𝑉1𝑘 ∙ 𝐴1 ∙ sin(2𝜋𝜔1 𝑡 + 𝛼1 )
(11)
⃗⃗⃗⃗⃗⃗
𝑉1𝑘 は k 目の骨格の 1 番目のモードの変位を表す。ここで振動数𝜔1 = 1、位相𝛼1 = 0
とし、これを（11）式に反映させると、
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
𝑊𝑘 (𝑡) = ⃗⃗⃗⃗⃗⃗
𝑉1𝑘 ∙ 𝐴1 ∙ sin(2𝜋𝑡)
1
3
となる。
（12）式より 2 𝜋, 2 𝜋, … ,
1
2
(12)
+ 𝑛 𝜋, 𝑛 ∈ 𝑵となる t において、骨格は最大の変位
量をとる。今回は 1 フレームごとの刻み t＝0.01 としており、25 フレーム目の t＝
0.25 のときが最大変位となるので、このときの変位を Blender 上で取得した。具体
的には、ヘモグロビンの原子を球のオブジェクトとして表示させたモデルに骨格を
設定し、25 フレーム目における原子の変位を得た。今回は振幅𝐴1 を 1,5,10 として、
骨格の動きを決定し、すべての変位について求めた。求められた変位と基準振動解
析結果の変位の誤差について、平均と RMSD を計算し、4.1.2 節の 表 2 にまとめ、
図 11 に示した。また、今回利用した最低振動数 10 個のモードの誤差についても、
平均と RMSD を計算し、図 12 に表した。同様に、二次構造法の骨格を原子のモデ
ルに適用させて、誤差を計算し、4.3 節の図 15 に表した。
24
4. 結果・考察
3 章で説明した提案手法によって作成された骨格を 3DCG ソフトウェアである
Blender22を用いて、タンパク質モデルに設定した結果を説明する。
4.1 サブドメイン法の骨格による結果
4.1.1
骨格の定性的評価
Blender 上で、サブドメイン法によって決定した骨格をタンパク質分子のモデル
（ヘモグロビン）に設定した結果を 図 9 に示す。図 9 において、灰色の楔形が骨格
である。モデルが設定した骨格によって動き、タンパク質の熱揺らぎが表現されて
いることが確認された。
図 9. サブドメイン法の骨格をヘモグロビンのモデルに設定した結果
gray wedge：骨格、yellow ribbon：リボンモデル、red ball & stick：ヘム鉄
22 http://blender.jp/
25
サブドメイン法の骨格によって動くタンパク質のモデルにおいて、構造生物学的
には安定とされる二次構造の中心部分が伸縮するような動きが観察された。これは
サブドメイン法による骨格の動きは平行移動のみで表されていて、タンパク質のへ
リックスやシートのような二次構造の中心付近を支点として傾くような振動の場
合、この二次構造を動かすのに少なくとも 2 つの骨格を必要としていたからと考え
られる（図 10）。今回は分割数が少なかったためこのような問題が生じたが、分割数
を増やし骨格の数を増やすことで、必然的にモデルの二次構造を動かす骨格が増え、
解消すると予想される。しかし、先にも述べたが構造的には安定な二次構造を 3DCG
ソフトで表す際には、二次構造に対して骨格を 1 つ決定する二次構造法のほうが直
感的には理解を得やすいと考える。
図 10. 骨格によるモデルの二次構造変形の仮説
26
4.1.2
骨格の定量的評価
3.3 節の方法で求められた変位と基準振動解析結果の変位の誤差を計算し、表 2 に
まとめ、図 11 に各値をプロットした結果を示す。Blender において、骨格による分
子モデルの動きが正確かどうかを確かめるために、振幅を 1,5,10 倍として誤差を計
算した。基準振動と骨格による動きの原子の誤差は線形に増加していることから、
骨格に動きを与える際の計算が Blender におけるフレーム間の補間が正しく行われ
ていたと考えられる。また、タンパク質の原子モデルに動きを設定する際に要した
時間の平均は、基準振動解析結果を反映させる場合約 890 秒、サブドメイン法の骨
格の場合約 12 秒であった。したがって、約 74 倍動きの設定を短縮させることがで
き、骨格によって軽処理化が実現できたと言える。
表 2 より炭素－水素間の結合距離が約 1 Å であるので振幅を解析結果の 10 倍に
して表示しても、化学的にはそれ以下の差である（平均 0.334 Å, RMSD 0.383Å）
。
各モードの誤差（図 12）は、平均で約 0.1 Å であることがわかった。骨格は全ての
振動モードを考慮して決定したので、特定の振動モードにおいては変位の差が大き
くなる可能性があった。しかし、その誤差におおきな差がなかったことから、1 つの
骨格で全ての振動モードの変位を十分に表現できていると考えられえる。
Blender において表現するときは、視覚的に振動をわかりやすくするため、振幅
を 5 倍程度にする。各モードの誤差の平均が約 0.1 Å であるが（図 12）、振幅を 5 倍
にしてもその差は 0.5 Å 程度であり、Blender において 0.5 Å は視覚的に微小な差で
あり、解析結果を十分に再現できていると言える。
表 2．ヘモグロビンの最低振動数 1 番目のモードにおける骨格による原子の変位
と基準振動解析結果との変位の誤差
計算式
振幅 1 倍
振幅 5 倍
振幅 10 倍
0.033Å
0.166 Å
0.334 Å
0.038 Å
0.191 Å
0.383 Å
∑𝑛𝑖=1 ∆𝑒1𝑖 (𝑡)
𝑛
(平均)
∑𝑛 ∆𝑒 (𝑡)
√ 𝑖=1 1𝑖
𝑛
(RMSD)
27
平均
RMSD
0.45
0.4
Å(=1.0x10-10m)
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0
2
4
6
8
振幅の倍率
図 11. 振幅の倍率による誤差の変化
28
10
12
RMSD
平均
0.14
Å(=1.0x10-10m)
0.12
0.1
0.08
0.06
0.04
0.02
0
0
2
4
6
8
最低振動数モードの番号
図 12. 最低振動数モードごとの誤差
29
10
12
4.2 二次構造法の骨格による結果
4.2.1
骨格の定性的評価
二次構造法で作成した骨格をヘモグロビンのモデルに設定した結果は図 13（A）
である。また、今回はへリックスが豊富なヘモグロビンとは対照的に、シートや
二次構造でない領域が多いアミラーゼ（PDB ID: 1AMY）についても骨格を作成し
モデルに設定した（ 図 14 ）
。
図 13（B）において、サブドメイン法の骨格で変形していたモデルの二次構造
（へリックス）は、二次構造法による骨格において変形は確認されなかった（図
13 A）。また二次構造法の骨格には平行移動だけでなく回転の動きも与えている
が、ここでも大きく変形した領域は確認されなかったことから、問題なく回転が計
算されていると考えられる。
B
A
図 13．設定された骨格の比較
（A）二次構造法による骨格をモデルに設定した結果。gray wedge；骨格、
yellow ribbon；リボンモデル、red ball & stick；ヘム鉄。
（B）サブドメイン法
による骨格をモデルに設定した結果。gray wedge；骨格、yellow ribbon；リボ
ンモデル、red ball & stick；ヘム鉄。
30
二次構造が少ないタンパク質であるアミラーゼについて二次構造法の骨格は、十
分に決定されていた（ 図 14 ）
。また骨格の動きによって、歪に変形する領域も確認さ
れなかった。そして本来は大域的な領域の変動しか観察されない基準振動を合成す
ることにより、局所的な変化（図 14、右下のループ領域）を確認することができた。
このようなループ領域の変化は、タンパク質の機能について考察することに役立
つ。アミラーゼが基質である多糖を分解後に、このループ領域によって送り出して
いるという仮説を立てることができる。
図 14. 二次構造法の骨格をアミラーゼのモデルに設定した結果
blue ribbon；アミラーゼのリボンモデル、gray wedge；骨格
31
4.3 サブドメイン法と二次構造法の比較
サブドメイン法と二次構造法の骨格をヘモグロビンの原子モデルに設定した際の、各
原子の変位と、基準振動解析の変位との差の平均を各モードについて計算した結果が 図
15 である。また、誤差に応じてヘモグロビンのリボンモデルを色分けし、図 16 と図
17 にて示した。図 16 がサブドメイン法、図 17 が二次構造法について、それぞれの
誤差を色分けした図になる。ここでは、骨格数は 2 つの方法で統一し、56 個とした。
図 16 ではほとんどの原子の誤差が 0.1 Å 以下であるのに対し、図 17 では 0.1 Å よ
り大きくなっていることが確認された。また、両手法とも一部のへリックスにおい
ては、誤差が 0.2 Å よりも大きくなることが確認された。図 15 より基準振動解析結
果との誤差は、二次構造法の方が大きかったが、二次構造法の各モードの誤差を平
均すると約 0.16Å であり、Blender において振動を見やすくするため、振幅を 5 倍に
した骨格によって原子を動かしても、その誤差は約 0.8 Å となることが予想された。
この差は、タンパク質全体のアニメーション表示において、視覚的に変化を確認で
きない程度であるが、原子モデルを表示した際には化学結合が不自然になる可能性
があるので注意すべきである。誤差の原因としては、二次構造法は時間刻み t によ
って得られる 3 点の変位をもとに、骨格の位置と角度を計算しているが、変化量が
微小すぎて角度の計算がうまく行われない部位のあることが考えられる。
また、現在は分子内の原子間の距離や衝突を考慮していないので、振幅の倍率原
子によっては化学的に無理な配座を取っている可能性する必要がある。
今回は 2 つの手法で骨格を作成したが、適用するモデルの表示方法によってその
優位性は異なると考える。タンパク質の二次構造をリボンで表したリボンモデルの
場合は、構造生物学的には安定とされる二次構造付近に骨格が決定されているほう
が感覚的にわかりやすく、二次構造法が良い。一方で、分子を構成する原子の電子
表面を表したサーフェスモデルでは、視覚的には二次構造を判断できないので、サ
ブドメイン法の骨格でも問題ないと考えられた。タンパク質のモデルについて骨格
を作成する場合は二次構造法のほうが一般的な生物学的な議論に向いているが、計
算が簡単であるサブドメイン法は汎用性が高く、二次構造が崩れている場合に有利
と思われる。
近隣で同方向に動く C原子をクラスタリングすることによって見出される原
子の集団運動は、タンパク質の二次構造と関連性があると言われている[14]。一方で、
二次構造によっては、必ずしもまとまった集団運動を示さない場合があることも言
32
われている。また、タンパク質の活性中心のような機能的に重要な部位は、原子の集
団運動に伴う二次構造の変化にあまり影響しないということが示唆されている。こ
れは、活性部位のような重要部位が、構造の動きによって変化してしまうと酵素と
しての機能が失われてしまうと考えると、妥当であるとされている。
このようにして、タンパク質を構成する原子の動きと立体的な近さをもとにした
クラスタリングは、原子の集団運動を特定し、構造的な特徴や機能を知る際 に利用
される。
サブドメイン法の平均
0.2
二次構造法の平均
0.18
0.16
Å(1.0x10-10)
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
0
2
4
6
8
10
最低振動数のモード番号
図 15. サブドメイン法と二次構造法の骨格による誤差の比較
33
12
B
A
C
D
■𝑒2𝑖 ≤ 0.05 Å ■0.05 Å < 𝑒2𝑖 ≤ 0.1 Å ■0.1 Å < 𝑒2𝑖 ≤ 0.2 Å ■0.2 Å < 𝑒2𝑖
図 16. サブドメイン法の骨格（第 2 振動モード）による誤差の色分け
𝑒2𝑖 ; 骨格の振幅が最大となるフレームにおける、最低振動数 2 番目のモードの i
番目の原子の誤差。（A）正面図。（B）
（A）から見たときの右側面図。（C）
（A）
から見たときの左側面図。（D）（A）から見たときの背面図。
34
B
A
C
D
■𝑒2𝑖 ≤ 0.05 Å ■0.05 Å < 𝑒2𝑖 ≤ 0.1 Å ■0.1 Å < 𝑒2𝑖 ≤ 0.2 Å ■0.2 Å < 𝑒2𝑖
図 17. 二次構造法の骨格（第 2 振動モード）による誤差の色分け
𝑒2𝑖 ; 骨格の振幅が最大となるフレームにおける、最低振動数 2 番目のモードの i
番目の原子の誤差。（A）正面図。（B）
（A）から見たときの右側面図。（C）
（A）
から見たときの左側面図。（D）（A）から見たときの背面図。
35
4.4 作成したアニメーションの生物学的な考察
分子アニメーションは、タンパク質の構造変化による機能の仮説を示すために利
用されてきた。タンパク質の構造変化や機能についての研究は、数多くの研究者に
よってなされ少しずつ明らかになってきてはいるが、唯一無二の構造変化が提唱さ
れているわけではない。そのため、多くの仮説が存在し、それらの仮説を表現する
ためにアニメーションが利用される。ここでは、本研究で取り扱ったヘモグロビン
およびアミラーゼについて、それらの機能の仮説がアニメーションで表現されてい
るという例を挙げる。そして、本研究で作成したアニメーションがどのように生物
学的に考察されるのかを述べる。
ヘモグロビンは血中の赤血球に含まれるタンパク質で、酸素分子と結合し、肺か
ら全身へと酸素を運搬する役割を担うタンパク質である。ヘモグロビンはサブユ
ニットとサブユニットと呼ばれる、2 種類のサブユニットそれぞれ 2 つから成る 4
量体のタンパク質である。また、ヘモグロビンは酸素を結合していない T（taut）
状態と結合した R（relax）状態が存在する。そして、それぞれのサブユニットが
酸素を結合し T 状態から R 状態に変化するわけではなく、1 つのサブユニットが酸
素を結合すると残りのサブユニットも R 状態へと変化し、酸素を結合しやすい構造
を取るという仮説が提唱され、アニメーションで表現されている 23。本研究で作成
したヘモグロビンの熱揺らぎを表現するアニメーションにおいては、サブユニッ
トとサブユニットそれぞれが対称的に、全てのサブユニットの構造が変化する動
きが観察された。これは、T 状態から R 状態へ全体が構造変化するという仮説に近
い動きであった。
アミラーゼは澱粉中に含まれるアミロースやアミロペクチンに代表される多糖を
分解するタンパク質である。このアミラーゼは工業的に利用価値が高く、水飴やブ
ドウ糖製造時の澱粉液化に利用されている 24。特にアミラーゼは糖鎖のグリコシド
結合を加水分解する酵素の総称で、作用する部位によってアミラーゼ、アミラー
ゼ、グルコアミラーゼに分けられる。糖鎖を用いた分解のアニメーションも存在
し、異なる分解作用の理解の助けとなっている 25。本研究では、アミラーゼの図
14 の右下のループ領域が大きく振れていることが、アニメーションによって観察
23 http://biochem.web.utah.edu/iwasa/projects/he moglobin.html
24 http://www.hbi-enzymes.com/HBI_Amylase.htm
25 http://www.youtube.com/watch?v=HTJy3lnPmL8
36
された。アミラーゼは多糖を分解する酵素であり、観察された大きく動くループ
領域が多糖を分解後に放出しているのではないかと推測された。このような、タン
パク質の熱揺らぎによる瞬間的な構造変化は、タンパク質の機能に関連するのでは
ないかという仮説が提唱されている[15]。
熱揺らぎによる構造変化を眺めることで、これまで考察されてこなかった、一時
的に顕れる特徴的な構造変化を捉えることができ、新たな実験計画の立案につなが
る。熱揺らぎのアニメーションは新たな仮説立案を支援し、タンパク質の機能理解
に役立つ可能性がある。
37
4.5 適用の拡大に関する考察
今回は k-平均法 を原子の領域決定に用いたが、分割数𝐾を決めるのは難しい問題であ
る。今後多くのタンパク質について骨格を決定するたび、クラスタ数𝐾を任意に決定す
る方法では時間を要する。そこで、分割数𝐾を決める指標について考察する。3DCG ソ
フトウェアにおいて、ボーンの数が少ないとモデルが取り残されるという問題が生じ
る。サブドメイン法、ならびに二次構造法では k-平均法によって原子を 𝐾個のサブド
メインに分割しているが、1 つの指標として、取り残されるポリゴンが無いように、
骨格数𝐾を設定する必要がある。タンパク質にはアミノ酸配列の類似性に応じて分類
される、ファミリーと呼ばれるものがある。例えば、これらのファミリーに代表され
るタンパク質に応じて、取り残されるポリゴンが無くなるときの骨格数を調べる。
そして、おおよそ何個の骨格が必要かを決めることができれば、同じファミリーに
属すタンパク質についても同様か近い分割数が妥当になると考えられる。このよう
に、タンパク質の分類、原子数、構造の特徴などを指標にして、分割数を予想できる
可能性がある。k-平均法において、情報量を最大にする指標を用いて分割数𝐾を決め
る手法が x-means 法として提案されている[16]。サブドメイン法において適用を試
みたが、ヘモグロビンで 700 程度と多くなる傾向があり、本研究においては十分な
実用性を得られなかった。
Blender を用いて、ヘモグロビンを構成する全原子 1 つ 1 つを球のオブジェクト
として取り込み、振動モードの 1 つを表示させようとしたところ、処理が完了する
まで数分から数十分の時間を要した。対して、分子ビューアーにて PDB ファイルか
ら読み込んだヘモグロビンの構造をリボンモデルや分子表面モデルに変換、ポリゴ
ンファイルを作成し Blender にて取り込んだ場合は数秒で読み込みが終了する。そ
れはタンパク質によっても異なるとはいえ、リボンモデル、分子表面モデルにおい
ても数個程度のオブジェクトで全体表示が可能であるからと考えられる。これは、
骨格の情報ファイルである BVH ファイルを取り込み、分子のモデルに設定した場
合も同様で、全原子の振動モードを表示させるより早い。また Blender でのモデル
を操作する場合でも、表示するオブジェクトが少ないほど軽処理で操作できる。も
し原子モデルを動かしたい場合は、本章でこれまで述べてきたように、タンパク質
の主鎖の一部の原子を表示させ、それに骨格を設定すれば良い。
38
4.6 今後の展望
本研究ではタンパク質の基準振動解析結果を用いて骨格を作り、これを用いて
3DCG ソフトウェアで利用されるモデルを動作させた。タンパク質だけではなく、
他の生体内高分子である糖や核酸の基準振動解析の結果があれば、本提案手法を利
用したアニメーションが作成できると言える。糖や核酸には、タンパク質のように
二次構造は無いが、特徴的な骨格を定義して二次構造法を適用するか、構成する原
子をクラスタリングによって分割し、骨格を作成するサブドメイン法であれば適用
できると考えられる。また、本研究では基準振動解析の結果を用いたが、時系列ごと
に原子の変位が記述された分子動力学シミュレーションのトラジェクトリデータを
利用することも可能である。まず、基準振動解析のデータによって骨格を作成して
おき、骨格を構成する原子についてトラジェクトリ座標から骨格の座標を再計算す
れば、トラジェクトリに沿った骨格の動きが得られる。すなわち、モーションキャプ
チャが容易に実現できる。基準振動解析が無い場合であっても、トラジェクトリを
主成分分析する Essential Dynamics の手法[17]で、基準振動解析と同様の振動モー
ドを取得することにより、同様にして本手法に用いることができると考えられる。
39
5. おわりに
これまで分子アニメーション制作で利用されるモデルの動きは、科学的根拠に基
づかない手作業や大規模計算機による計算によって与えられることが多かった。本
研究ではそれらの問題点に対して、分子の熱揺らぎを計算する基準振動解析に基づ
いた骨格を作成して用いる手法を提案した。この手法を用いてモデルを軽処理かつ
生体内で予想される動きで動かすことができた。制作を行うにあたって、骨格ファ
イルを扱いやすいように何種類かのタンパク質について骨格ファイルを作成しデー
タベース化することで利用が普及すると考えられる。
本格的な分子アニメーションによる可視化の敷居が下がってきている。しかし、
分子の立体構造や動きを評価するのは実験生物学者であり、総合的な検証を繰り返
す必要がある。特にタンパク質間の相互作用の可視化は、原子間の分子間力などの
多様なパラメータを考慮した上で、現実的な動きが再現できていなければいけない。
また科学的根拠のあるシミュレーション結果を利用する分子アニメーションも、必
ずしも正しいとは限らない。したがって、評価は個別の研究によってそれぞれ評価
されるべきであり、情報学的な可視化手法はそういった評価の積み重ねを実績とし
て、信頼性を高めていかなければならない。そのためには、今回のような代表的なタ
ンパク質での成果のみならず、今後多くの種類のタンパク質においても同様の検証
をしていく必要がある。
また、計算機の処理能力と 3DCG ソフトウェアの発展を考えると、計算機がより
複雑な計算を短時間に処理する能力を獲得することは容易に考えられる。計算機の
中央演算処理装置である CPU の開発は、近年の半導体の高密度化と微細化によって
数 nm プロセスという回路における物理的な限界に近づきつつあり、省電力や小型
化といった他の特徴に注力しつつある。一方で、特に GPU などのグラフィック分野
における処理能力はいまだに増加傾向にあり、ゲームや映像分野の発展も相まって
高処理能力化すると考えられる。しかしながら、計算量が多い 1 原子単位での分子
動力学計算や熱揺らぎの動きを 3DCG ソフトウェアで容易に再現するには、年月を
要すと考えられる。そうしたときに、本研究で提案した骨格を用いた手法のような
計算量を軽量化する方法は、未来においても陳腐化せず、必要とされると想定され
る。そのためには、必要とされる理想的な分子アニメーションを常に意識し、それを
現状の環境で実現するための簡易化モデルを考えていかなえればならない。
40
謝辞
本研究の機会を与えてくださり、研究に必要な技術や理論に関する献身的な御指
導と多大なる御助言、御援助を賜りました、独立行政法人産業技術総合研究所 健
康工学部門 上野豊 主任研究員に心より厚く御礼を申し上げます。
研究を進めるにあたり、丁寧な御指導、御鞭撻を頂戴致しました、主指導教官で
ある情報科学研究科 金谷重彦 教授に深く感謝の意を表します。
本研究を遂行するにあたり、有益な御助言を頂戴致しました、情報科学研究科
佐藤嘉伸 教授に厚く御礼申し上げます。
本研究に関して様々な御意見を頂きました、情報科学研究科 Md. Altaf-ULAmin 准教授に心より感謝の意を表します。
本研究を取り組むにあたり、日頃より多岐にわたる御助言、御指導を賜りまし
た、情報科学研究科 佐藤哲大 助教に深く感謝致します。
本研究に関して様々な御意見を頂きました、情報科学研究科 小野直亮 助教に心
より感謝を申し上げます。
本研究の遂行にあたり、多くの御助言を頂戴致しました、情報科学研究科 杉浦
忠男 准教授に深く感謝致します。
本研究に関して様々な御意見を頂きました、情報科学研究科 西岡孝明 特任教授
に心より感謝を申し上げます。
本研究を遂行するにあたり、様々なご意見を頂きました、情報科学研究科 黄銘
特任助教に深く感謝いたします。
本研究を進めるにあたり多くの御援助を頂きました独立行政法人産業技術総合研
究所 健康工学部門の皆様に深く感謝致します。
本研究に関して様々な御意見と御助言を頂きました、前情報科学研究科 池田俊
様、桂樹哲雄 様と計算システムズ生物学研究室MRグループの皆様に心より感謝致
します。
そして、充実した学生生活を送ることに大きく貢献して下さった計算システムズ
生物学研究室の皆様に心より感謝致します。
41
参考文献
[1] B. Alberts, A. Johnson, J. Lewis, D. Morgan, M. Raff, K. Roberts, and P. Walter,
“Molecular biology of the cell,” 5th ed., Garland Science, (2007).
[2]
須藤健太郎, 越塚誠一, “粒子法による泡の運動シミュレーション,” 第 19 回計
算力学公演論文集, pp.703-704 (2006).
[3]
G. T. Johnson, L. Autin, D. S. Goodsell, M. F. Sanner, and A. J. Olson, "ePMV
Embeds Molecular Modeling into Professional Animation Software Environments,"
Structure. Vol. 19, Issue 3, pp.293-303 (2011).
[4]
Y. Porozov, R. Andrei, and M. Zoppe, “Visualization of moving biomolecules: a
new approach based on professional,” Network Tools and Applications in Biology,
pp12-15 (2007).
[5]
T. A. Moller, E. Haines, and N. Hoffman, “Real-Time Rendering,” 3 rd ed., A K
Peters, Ltd., (2008).
[6]
M. Maddock, and S. Maddock, “Motion capture file formats explained,” Department of
Computer Science, University of Sheffield, (2001).
[7]
M. F. Sanner, “Python: a programming language for software integration and
development,” Journal of Molecular Graphics and Modeling. Vol. 17, pp.57-61 (1999).
[8]
古明池勇人, 上林正巳, 長嶋雲兵, “生体分子の分子動力学シミュレーション,”
Journal of Chemical Software. Vol. 16, No. 1, pp.1-36 (2000).
[9]
中村周吾, “核酸・タンパク質複合体の基準振動モード解析”, 計測と制御, 第 40
巻, 第 8 号, pp.598-603 (2001).
[10]
H. Wako, and S. Endo, “Normal mode analysis based on an elastic network model for
biomolecules in the Protein Data Bank, which uses dihedral angles as independent
variables,” Computational biology and Chemistry. Vol. 44, pp.22-30 (2013).
[11]
H. Wako, M. Kato, and S. Endo, “ProMode: a database of normal mode analyses on
protein molecules with a full-atom model,” Bioinformatics. Vol. 20, No. 13, pp.20352043 (2004).
[12]
平井有三, “はじめてのパターン認識”, 森北出版株式会社, pp.155-157 (2012).
42
[13]
W. Kabsch, and C. Sander, “Dictionary of protein secondary structure: pattern
recognition of hydrogen-bonded and geometrical features,” Biopolymers. Vol. 22, No.
12, pp.2577-2637 (1983).
[14]
鎌田真由美, 戸田幹人, “タンパク質分子の構造ダイナミクス: ウェーブレット
変換による解析,” 数理統計 第 62 巻, 第 2 号, pp.203-220 (2014).
[15]
M. Kokkindis, N. M. Glykos, V. E. Fadouloglou, “Protein flexibility and enzymatic
catalysis,” Advances in Protein Chemistry and Structural Biology. Vol. 87, pp.181-218
(2012).
[16]
T. Ishikawa, ”Extended k-means with an Efficient Estimation of the Number of the
Number of Clusters,” Lecture Notes in Computer Science. Vol. 1983, pp.17-22 (2000).
[17]
S. Hayward, and B. L. de Groot, “Normal modes and essential dynamics,” Methods
in Molecular Biology. Vol. 443, pp.89–106 (2008).
43