特開 2022-102389 生体高分子の動的相関量解析方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2022102389

(43)【公開日】2022-07-07

(54)【発明の名称】生体高分子の動的相関量解析方法

(51)【国際特許分類】

   G01N 33/15 20060101AFI20220630BHJP

   G16C 10/00 20190101ALI20220630BHJP

【ＦＩ】

G01N33/15 Z

G16C10/00

【審査請求】未請求

【請求項の数】3

【出願形態】ＯＬ

(21)【出願番号】P 2020217087

(22)【出願日】2020-12-25

(71)【出願人】

【識別番号】000125347

【氏名又は名称】学校法人近畿大学

(74)【代理人】

【識別番号】100118924

【弁理士】

【氏名又は名称】廣幸 正樹

(72)【発明者】

【氏名】米澤 康滋

(57)【要約】      （修正有）

【課題】ＭＤでは異なる時刻の分子座標決定は従来法では誤差が不可避であり分子内原子群の動的な相関を精密に評価する事は不可能であった。
【解決手段】本発明は、ＭＤによって平衡状態になった生体高分子モデルの部分構造同士の動的相関を解析する分子動的相関量解析方法であって、前記平衡状態から連続するＮ個の計算時刻における前記部分構造を構成する原子の座標を得る工程と、前記計算時刻毎に個々の前記部分構造内の原子間距離を代表する主要距離を求める工程と、前記個々の部分構造の前記主要距離を確率分布に変換する工程と、前記部分構造の全ての２個の組み合わせに対して、それぞれの前記部分構造の主要距離から同時確率分布に変換する工程と、前記それぞれの部分構造の確率分布と前記同時確率分布から、前記それぞれの部分構造同士の動的相関として相互情報量を求める工程を有することを特徴とする分子動的相関量解析方法である。
【選択図】図２

【特許請求の範囲】

【請求項1】

生体分子の分子動力学シミュレーション（ＭＤ）によって平衡状態になった生体高分子モデルの各部分構造同士の動的相関を、コンピュータを用いて解析する生体高分子の動的相関量解析方法であって、
前記平衡状態から連続するＮ個の計算時刻における前記部分構造を構成する原子の座標を得る工程と、
前記計算時刻毎に個々の前記部分構造内の原子間距離を代表する主要距離を求める工程と、
前記個々の部分構造の前記主要距離を確率分布に変換する工程と、
前記部分構造の全ての２個の組み合わせに対して、それぞれの前記部分構造の主要距離から同時確率分布に変換する工程と、
前記それぞれの部分構造の確率分布と前記同時確率分布から、前記それぞれの部分構造同士の動的相関として相互情報量を求める工程を有することを特徴とする生体高分子の動的相関量解析方法。

【請求項2】

前記主要距離を求める工程は、
前記部分構造を構成する原子間距離を求める工程と、
前記Ｎ個の計算時刻にわたる前記原子間距離の平均となる時間平均距離を求める工程と、
前記原子間距離と前記時間平均距離の差である偏差を求める工程と、
前記部分構造内の原子間距離の２つの組み合わせに対して共分散を求める工程と、
前記部分構造内の全ての共分散を要素とする共分散行列を求める工程と、
前記共分散行列の固有ベクトルの内最大の固有ベクトルを前記主要距離とする工程であることを特徴とする請求項１に記載された生体高分子の動的相関量解析方法。

【請求項3】

さらに、前記部分構造同士の前記相互情報量を２次元で表した相互情報量マップを作成する工程を含むことを特徴とする請求項１または２の何れかに記載された生体高分子の動的相関量解析方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、コンピューターを用いたシミュレーションにより生体高分子の動的相関量解析に関する。

【背景技術】

【0002】

タンパク質等の生体高分子内部では様々な原子が密接に関連して分子機能を制御している。近年この生体高分子内部の原子の動きを、シミュレーションを用いて把握することで、創薬やタンパク質合成の開発が進められている。

【0003】

特許文献１には、標的タンパク質と薬物分子との結合過程における分子運動を解析する方法が開示されている。ここでは、（１）標的タンパク質と薬物分子との複合体の立体構造について回転対称境界条件下における平衡化を行う工程、（２）工程（１）で得られた平衡化構造について分子動力学シミュレーションを行う工程、および、（３）工程（２）で得られた構造と初期構造との根平均２乗偏差を配座成分、並進成分、回転成分へ分解する工程、を含む方法により分子運動を解析する。

【先行技術文献】

【特許文献】

【0004】

【特許文献1】特開２００２－２４３７２２号公報

【発明の概要】

【発明が解決しようとする課題】

【0005】

タンパク質等、生体分子の分子動力学シミュレーション（以後「ＭＤ」と呼ぶ。）は、原子の位置の精密な情報を含み、動的に変化する分子構造の時系列データを生成する。これまでＭＤのデータ解析は主に特許文献１のように最小二乗法を使った分子重ね合わせ方法を使ってなされてきた。

【0006】

生体高分子は熱揺動によってその構造を柔軟に時々刻々変化させるので分子内原子の動的な絶対座標を精密決定する事は原理的に不可能である。その為、これまで分子内の局所部分の座標値は重ね合わせによる平均位置として決定された。

【0007】

しかし、異なる時刻の分子座標決定は誤差を避ける事が出来ず分子内原子群の動的な相関関係を精密に評価する事は不可能であった。そのため、活性中心以外の部分の活動が酵素活性や触媒活性を増減するアリステリック効果をシミュレーションで検討することはできなかった。

【課題を解決するための手段】

【0008】

本発明は上記の課題に鑑みて想到されたものであり、ＭＤで平衡状態になった計算結果において、原子間の距離のみの変動を調べることで、任意の生体高分子を構成する部分構造間（タンパク質を構成するアミノ酸同士やＤＮＡの塩基同士）の動的相関を解析するものである。

【0009】

より具体的に本発明に係る生体高分子の動的相関量解析方法は、
生体分子の分子動力学シミュレーション（ＭＤ）によって平衡状態になった生体高分子モデルの各部分構造同士の動的相関を、コンピュータを用いて解析する新規解析方法であって、
前記平衡状態から連続するＮ個の計算時刻における前記部分構造を構成する原子の座標を得る工程と、
前記計算時刻毎に個々の前記部分構造内の原子間距離を代表する主要距離を求める工程と、
前記個々の部分構造の前記主要距離を確率分布に変換する工程と、
前記部分構造の全ての２個の組み合わせに対して、それぞれの前記部分構造の前記主要距離から同時確率分布に変換する工程と、
前記それぞれの部分構造の確率分布と前記同時確率分布から、前記それぞれの部分構造同士の動的相関の指標として相互情報量を求める工程を有することを特徴とする。

【発明の効果】

【0010】

本発明に係る生体高分子の動的相関量解析方法によれば、生体高分子を構成する部分構造内の原子間距離だけを用いるため、計算時刻間の生体高分子モデル同士を、重ね合わせを用いず動的な解析を行うことができるという効果を奏する。その結果、タンパク質同士が接近、結合した際のアリステリック効果を評価することができ、新たな薬剤の開発に寄与することができる。

【0011】

また、従来用いていた分子動力学シミュレーションのソフト資産をそのまま利用することができるので、導入もしやすいという効果を奏する。

【図面の簡単な説明】

【0012】

【図1】本発明に係る生体高分子の動的相関量解析を行うハード構成を示す図である。

【図2】本発明に係る生体高分子の動的相関量解析の概略を示す図である。

【図3】本発明に係る生体高分子の動的相関量解析を行うフローチャートである。

【図4】本発明に係る生体高分子の動的相関量解析を行うフローチャートである。

【図5】生体高分子の動的相関量解析の詳細な手順を示す図である。

【図6】アミノ酸中の原子から座標、原子間距離、時間平均距離を説明する図である。

【図7】偏差、偏差の共分散、共分散行列、距離固有ベクトル、主要距離を説明する図である。

【図8】主要距離を説明する図である。

【図9】確率分布および同時確率分布を説明する図である。

【図10】ＰＤＺのリボン図である。

【図11】ＰＤＺにリガンドペプチドが結合されていない場合と、結合している場合のＰＤＺの各アミノ酸同士の動的相関量を実際に求めた出力図（相互情報量マップ）である。

【発明を実施するための形態】

【0013】

以下に本発明に係る生体高分子の動的相関量解析方法について図面および実施例を示し説明を行う。なお、以下の説明は、本発明の一実施形態および一実施例を例示するものであり、本発明が以下の説明に限定されるものではない。以下の説明は本発明の趣旨を逸脱しない範囲で改変することができる。

【0014】

＜ハード構成＞
本発明に係る生体高分子の動的相関量解析方法は、コンピュータを用いて行われる。図１には、ハード的な構成を示す。コンピュータ１は中央処理装置（ＣＰＵ：Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｅｒ Ｕｎｉｔ）１０とメモリ１２で構成され、周辺機器として、入出力装置が備えられる。入出力装置は、一般に、キーボード、マウスといった入力デバイス１６と、ディスプレイ１８といった出力装置が利用できる。また、外部メモリ１４が備えられていてもよい。外部メモリ１４としては、ハードディスクドライブやソリッドステートドライブ、クラウドメモリ等が好適に利用できる。

【0015】

メモリ１２中には、適当なＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）と、アプリケーションソフトがロードされる。生体高分子の動的相関量解析方法では、ＭＤの出力を用いるので、アプリケーションとしてはＭＤと、生体高分子の動的相関量解析用のソフトウエア（ＹＺとする）が備えられる。ＭＤは利用可能なものが種々あるが、特に限定されるものではない。Ａｍｂｅｒ、Ｇｒｏｍａｃｓ等の世界的に用いられている主要プログラムを適宜利用することができる。

【0016】

解析に供される生体高分子は、タンパク質やＤＮＡ分子が対象となる。これらの構成する部分構造とは、タンパク質の場合は、アミノ酸であり、ＤＮＡ分子の場合は、塩基である。部分構造は、単一の物質（アミノ酸やＤＮＡ）でなくてもよく、連続する複数のアミノ酸やＤＮＡを１つの部分構造として扱ってもよい。また、複数のタンパク質やＤＮＡの塊がある場合は、それぞれのタンパク質やＤＮＡを部分構造としてもよい。また、タンパク質に化合物が結合したものであってもよい。解析対象となる生体高分子の原子座標に関する情報は、予め結晶法やＸ線回折、ＮＭＲ（核磁気共鳴）といった方法で取得されている。

【0017】

解析は、ＭＤによる生体高分子の平衡状態の算出と、平衡状態からＮ個の連続した生体高分子座標の取り出し、そして、生体高分子を構成する部分構造同士の動的解析という手順で進む。なお、以後の説明は生体高分子としてタンパク質を、部分構造としてアミノ酸を例として説明を行うが、生体高分子としてＤＮＡを、部分構造として塩基としても全く同じ手順で結果を得ることができる。

【0018】

よく知られているように、ＭＤは、物質を構成する原子同士の結合情報に基づいて、温度や原子同士に働く力に基づいてタンパク質の平衡状態を一定時刻毎に算出する。平衡状態を算出した時刻を「計算時刻」と呼ぶ。隣接する計算時刻の間の時間を「単位計算時間」と呼ぶ。単位計算時間は数ｆｓ（フェムト秒）毎程度が設定される。平衡状態とは、一定の時間の間の評価指数（例えばエネルギー）の変化が所定値以下になった状態で判断される。

【0019】

＜解析の流れ＞
本発明に係る生体高分子の動的相関量解析方法について、図２を用いて、その概略をまず説明する。解析対象となるタンパク質はＡ１からＡｎのｎ個のアミノ酸で構成されているとする。これらのアミノ酸を構成する原子の座標は与えられている。これをＭＤを用いてシミュレーションすることで、平衡状態を得ることができる。

【0020】

本発明に係る生体高分子の動的相関量解析方法は、ＭＤによるシミュレーションが平衡状態になった状態から得た連続したＮ個のモデルデータを使う（図２（ａ））。これは（単位計算時間×Ｎ）の間のモデルの時間変動のデータである。時刻をｔとすると、ｔ＝ｔ_１からｔ＝ｔ_Ｎまでのデータとなる。なお、モデルデータには、対象となるタンパク質を構成する原子の座標データが含まれる。

【0021】

本発明では、アミノ酸毎にこれら原子間距離を代表する「主要距離Ｌ」を定義する（図２（ａ））。主要距離Ｌは、その計算時刻の際のそのアミノ酸を構成する原子同士の原子間距離を代表できれば、限定はされない。例えば、アミノ酸内の全ての原子同士の和などである。なお、アミノ酸内の全ての原子同士の原子間距離に対して全ての計算時刻に渡る偏差を求め、距離同士の共分散を成分とする共分散行列の固有ベクトル成分の和を主要距離Ｌとすることもできる。このような主要距離Ｌの算出の一例の詳は後述する。

【0022】

各計算時刻におけるアミノ酸毎に主要距離Ｌを求めた後に、各アミノ酸毎に、Ｎ個の主要距離Ｌを確率分布に変換する（図２（ｂ））。これは主要距離Ｌを微小区間に分け、その区間内の密度を求める事で、度数に変換し、規格化することで求める事ができる。これは公知の方法である。主要距離Ｌに対する確率分布ｐは全てのアミノ酸について求められる。それぞれ、ｐ（Ａ１）、ｐ（Ａ２）、・・・、ｐ（Ａｎ）と表される。

【0023】

次に全てのアミノ酸から２つのアミノ酸を選び、それぞれのＮ個の主要距離を平面にプロットし、散乱分布を作り、同時確率分布に変換する（図２（ｃ））。それぞれのアミノ酸の主要距離Ｌに基づく確率分布は、同時確率分布の周辺確率分布となる。

【0024】

そして、これらの２つのアミノ酸の分子動的相関量を、それぞれの確率分布と同時確率分布から相互情報量Ｉｘｙを求めることで決定する（図２（ｄ））。相互情報量Ｉｘｙは、互いの確率分布の関連性を示す指標であり、互いに独立であればゼロとなり、互いの確率分布が同じであれば、エントロピーに一致する。図２（ｄ）では、選ばれたそれぞれの確率分布はｐ（Ａｘ）、ｐ（Ａｙ）で表され、同時確率分布はｐ（Ａｘ，Ａｙ）で示した。

【0025】

ここで求めた相互情報量Ｉｘｙは、Ｎ個の計算時刻における各アミノ酸内の原子間距離だけから求められたものであり、相互情報量Ｉｘｙが高ければ、それらのアミノ酸同士は、同期して動いているといってよい。モデルデータ中の全てのアミノ酸のペアに対して相互情報量Ｉｘｙを求め、平面にプロットすることで、解析対象となるタンパク質の分子動的相関量を求めることができる（図２（ｅ））。

【0026】

図２（ｅ）は、相互情報量Ｉｘｙをプロットしたもの（以下「相互情報量マップ」と呼ぶ。）を例示する。縦軸および横軸は、アミノ酸番号である。相互情報量Ｉｘｙが大きいものを黒丸の大きさで示した。この図では、Ａ２とＡｎで相互情報量Ｉｘｙが大きくなっており、これら１つのアミノ酸は共同して動いていることを示している。

【0027】

このように、本発明に係る生体高分子の動的相関量解析方法では、原子間距離を用いてアミノ酸同士の動きの相関性を求めることで、動的な相関量を解析する。以下詳細に説明する。なお、図３及び図４に本発明に係る生体高分子の動的相関量解析の処理フローを示す。以下図３及び図４のフローに従い各処理を詳説する。

【0028】

＜モデルデータの取得：ステップＳ１０２＞
図３の処理が開始されると（ステップＳ１００）、モデルデータが取得される（ステップＳ１０２）。図５にモデルデータの取得までの概略を示す。モデルとなるタンパク質には、Ａ１からＡｎまでｎ個のアミノ酸が連結しているとする（図５（ａ））。それぞれのアミノ酸には、順に番号が付されたｋ_１、ｋ_２、・・・、ｋ_ｎといったｋ_ｋ個の原子があるとする。つまり、アミノ酸Ａ１には、ｋ_１個の原子があり、アミノ酸Ａ２には、ｋ_２個の原子があり、アミノ酸Ａｎにはｋ_ｎ個の原子があるということである。

【0029】

また、アミノ酸番号を左上に記載し、［］をつけることで、原子の座標を表す。ａ１、ａ２、・・・、ａｋ_ｎでそれぞれの原子を表し、［ａｋ_１］、［ａｋ_２］、・・・、［ａｋ_ｎ］で、各原子の座標を表す。

【0030】

したがって、アミノ酸Ａ２には、^２ａｋ_１、^２ａｋ_２、・・・、^２ａｋ_ｎの原子があり、それぞれの原子の座標は［^２ａｋ_１］、［^２ａｋ_２］、・・・、［^２ａｋ_ｎ］と表される。例えば、［^３ａ５］ならば、Ａ３アミノ酸中の５番目の原子の座標を表す（図５（ｂ））。なお、アミノ酸内の原子は炭素、酸素、窒素といった大きな原子だけであってよく、水素は考慮しなくてもよい。

【0031】

この原子の座標データおよびその他の条件と共にＭＤで平衡状態をシミュレートされる（図５（ｃ））。図５（ｄ）には、シミュレーションの経過を示す。横軸は計算回数（計算時刻）であり、縦軸は収束の指標である。指標の変化が大きい部分は平衡化を行っている部分であり、指標の変化が落ち着くとサンプリング期間と呼ばれ、平衡化となった状態である。モデルデータは、このサンプリング期間からＮ個（計算時刻ｔ＝ｔ_１からｔ＝ｔ_Ｎまで）の連続する計算時刻の情報を取り出す（図５（ｅ））。

【0032】

隣接する時刻の間は単位計算時間であり、各計算時刻毎にｎ個のアミノ酸を有するモデルデータが存在する。アミノ酸番号の後ろに計算時刻を付す。例えば、Ａ２（ｔ_２）は、時刻ｔ_２の時のアミノ酸番号２のアミノ酸を示す。

【0033】

上述したように、アミノ酸Ａ２には、ｋ_２個の原子があり、それぞれの原子の座標をあらわに示すと、［^２ａｋ_１］（ｔ_２）、［^２ａｋ_２］（ｔ_２）、・・・、［^２ａｋ_ｎ］（ｔ_２）となる（図５（ｆ））。

【0034】

＜主要距離の算出：ステップＳ１０４＞
平衡状態からモデルデータを取得した後は主要距離Ｌの算出を行う（図３：ステップＳ１０４）。主要距離Ｌは、特定時刻の特定アミノ酸内の全ての原子間距離を代表できる数値であればよい。例えば、全ての原子間距離の和であってもよい。しかし、より妥当性の高い方法として、以下に詳細を述べる。処理のフローは図４のフローとなる。

【0035】

＜原子間距離：ステップＳ２００＞
主要距離Ｌの算出を説明するが、図６乃至図８も用いる。図６ではアミノ酸Ａｍには３つの原子（α、β、γ）があるとする。まず、原子間距離を定義する。計算時刻ｔにおけるｍ番目のアミノ酸Ａｍ内の任意の２つの原子間距離を^ｍｄｉｊ（ｔ）で表す。つまり、時刻ｔの時のアミノ酸Ａｍ中のｉ番目とｊ番目の原子間の距離は（１）式で表される。（１）式はｉ番目とｊ番目の原子間距離をユークリッド距離として求めたものであり、ルートの中は、単位軸成分の差の２乗の和を簡略的に記載した。

【0036】

【数1】

【0037】

図６（ａ）では、アミノ酸Ａｍ内に３つの原子があり、それぞれの座標が図６（ｂ）で表されることを示している。これはステップＳ１０２で得られたものである。なお、アミノ酸Ａｍの前にはアミノ酸Ａｍ－１がありｍの後にはアミノ酸Ａｍ＋１がある。アミノ酸の順は例えばＮ末端から数えた数としてよい。次に図６（ｃ）のように原子間距離^ｍｄｉｊ（ｔ）が求められる。これは（１）式に基づくものである。

【0038】

＜時間平均距離：ステップＳ２０２＞
次に、ｍ番目のアミノ酸のｉ番目とｊ番目の原子間の距離の時間平均を時間平均距離＜^ｍｄｉｊ＞とし、（２）式のように求める。図６では、図６（ｄ）のように、時間平均距離が求められる。これは（２）式に基づくものである。

【0039】

【数2】

【0040】

＜原子間距離の偏差：ステップＳ２０４＞
ｍ番目のアミノ酸の特定時刻ｔの原子間距離の偏差ＤＳ_ｉｊ（ｔ）を（３）式のように求める。これは、原子間距離と時間平均距離の差である。

【0041】

【数3】

【0042】

偏差ＤＳ_ｉｊ（ｔ）は、全ての距離について、また、全ての計算時刻に渡って求められる。図６（ａ）の例では、図７（ａ）のように求められる。図７（ａ）では、αβ間の距離、βγ間の距離、γα間の距離の３つ距離について、それぞれの偏差が計算時刻ｔ_１からｔ_Ｎまで求められる。それぞれの偏差はＤＳ_αβ、ＤＳ_βγ、ＤＳ_γαと示した。

【0043】

＜原子間距離の偏差の共分散：ステップＳ２０６＞
次にアミノ酸毎に、全ての原子間距離の偏差同士の共分散を（４）式のように求める。（４）式は、アミノ酸Ａｍ内において、計算時刻ｔの時の、ｉ番目とｊ番目の原子間距離の偏差と、ｋ番目とｌ番目の原子間距離の偏差の共分散を求め、全計算時刻について平均を取っている。この共分散を^ｍｄｉｓｔＰＣＡ_{ｉｊ／ｋｌ}とする。図７（ｂ）に同様の処理を示す。

【0044】

【数4】

【0045】

^ｍｄｉｓｔＰＣＡ_{ｉｊ／ｋｌ}は、ｍ番目のアミノ酸内の全ての原子間距離で全計算時刻ｔ_１からｔ_Ｎまでの偏差を用いて求められる。また、原子がｋ_ｋ個あるアミノ酸では、ｋ_ｋＣ２通りの原子間距離の偏差の組があり、その内の任意の２組の組み合わせについて原子間距離の偏差の共分散が求められる。この、^ｍｄｉｓｔＰＣＡ_{ｉｊ／ｋｌ}は、ｍ番目のアミノ酸の全計算時刻の原子間距離を代表する原子間距離に関する情報である。

【0046】

原子間距離の偏差の共分散は、全てのアミノ酸について、求められる。あらわに書くと、^１ｄｉｓｔＰＣＡ_{ｉｊ／ｋｌ}、・・・、^ｎｄｉｓｔＰＣＡ_{ｉｊ／ｋｌ}である。

【0047】

＜共分散行列：ステップＳ２０８＞
（４）式のように求めたｍ番目のアミノ酸の原子間距離の偏差の共分散の個々を要素とする共分散行列を（５）式のように作成する。（５）式はアミノ酸Ａｍの場合を例示しており、アミノ酸Ａｍについて、（５）式の共分散行列が１つ作成される。

【0048】

【数5】

【0049】

この共分散行列は、全てのアミノ酸について求められる。つまり、（５）式のような共分散行列は解析対象となるタンパク質に対してアミノ酸の個数分（ｎ個）作成される。

【0050】

図７（ｃ）では、原子が３つの場合（α、β、γ）のアミノ酸の共分散行列を表す。

【0051】

＜距離固有ベクトル：ステップＳ２１０＞
各アミノ酸に対して（５）式の共分散行列が作成されたら、この共分散行列の固有値と固有ベクトルを求める。例えば、（５）式は、ｋｍ×ｋｍの正方行列なので、固有値および固有ベクトルはｋｍ個求められる。これらの固有ベクトルを用いて、このアミノ酸における原子間距離の主成分軸を算出し、距離固有ベクトルとする。この距離固有ベクトルは、時刻ｔ_１からｔ_Ｎに渡って、該当するアミノ酸（例えば上記のＡｍ）内の全ての原子間距離の分散を反映している。

【0052】

主成分軸（距離固有ベクトル）の算出方法は、特に限定されない。例えば、固有ベクトルの内の最大の固有ベクトルを主成分軸としてもよい。しかし、共分散行列の性質をよく表す方法としては、固有値が大きな順に和を取り、その和が固有値の総和に対して一定の値以上になるように固有ベクトルを選び、それらを加えたものを主成分軸（距離固有ベクトル）とする方法が好適である。言い換えると、固有ベクトルの占有率が一定の値以上になるように固有ベクトルの和を取る。

【0053】

ここでの一定値は、５０％以上、９０％以下の範囲で決定するのが好適である。５０％未満では、全体を代表しているというには不十分であり、９０％より大きい値では、残り１０％を加えても主成分軸としては大きく変化しないからである。後述する実施例では、８０％以上になるように、固有ベクトルが選ばれ、主成分軸が求められる。

【0054】

＜原子間距離の主成分軸への射影：ステップＳ２１２、ステップＳ２１４＞
次に図８を参照する。図８は、アミノ酸Ａｍの計算時刻ｔ_１からｔ_Ｎまでを表す。各計算時刻において、原子は位置が異なり、得られる原子間距離は、計算時刻毎に変わる。ここで、（６）式に示すように、各アミノ酸の各計算時刻ｔにおいての原子間距離を主成分軸に投射し、その和をＰｒｏｊ_ｎ（ｔ）として求める。

【0055】

【数6】

【0056】

ここで、Ｅｉｇｅｎ＿ｖａｌｕｅ（ｎ，ｉ）は、ｎ番目のアミノ酸のｉ番目の固有値を表し、Ｅｉｇｅｎ＿ｖｅｃｔｏｒ（ｎ，ｉ）は、ｎ番目のアミノ酸のｉ番目の固有ベクトルを表す。Σの中は、固有ベクトルの占有率が一定値以上になるまで足し合わせることを示し、距離固有ベクトル（主成分軸）を示す。主成分軸は、１つのアミノ酸に対して、計算時刻にかかわらず１つ決まる。

【0057】

また、Ｄｎ（ｔ）は時刻ｔにおける、ｎ番目のアミノ酸の中で取りうる原子間距離を成分とするベクトルである。（６）式は、これらの内積を求めることで、主成分軸への射影を行っている。たとえば、Ｐｒｏｊ_２（ｔ_４）は、時刻ｔ４に於いて、アミノ酸Ａ２内の全ての原子間距離を主成分軸（距離固有ベクトル）へ投影し、その長さを合計したものである。

【0058】

したがって、１つのアミノ酸ｍに対してＰｒｏｊ_ｍ（ｔ_１）からＰｒｏｊ_ｍ（ｔ_Ｎ）までのＮ個（計算時刻の総数）のＰｒｏｊ_ｍ（ｔ）が求められる。このＰｒｏｊ_ｍ（ｔ）を各計算時刻の主要距離Ｌとする。全てのアミノ酸に対してＰｒｏｊ_ｎ（ｔ）を求めたら、処理はメインルーチンに戻される（ステップＳ２１６）。

【0059】

＜確率分布への変換：ステップＳ１０６＞
主要距離Ｌ（Ｐｒｏｊ_ｎ（ｔ））を度数分布と考えて、確率分布に変換する。これは公知の方法で行われてよい。例えば、一例を図９に示す。図９（ａ）を参照し、アミノ酸Ａｍの主要距離Ｌを横軸とし、横軸上に主要距離Ｌをプロットする。そして、主要距離Ｌの微小区間毎に密度を求め、縦軸にとる（図９（ｂ））。そして、密度を規格化することで確率密度に変換し、確率分布ｐ（Ａｍ）を得る（図９（ｃ））。確率分布ｐ（Ａｍ）は離散的であってもよい。この工程は図２（ｂ）に示したように全てのアミノ酸について行われる。すなわち、このような確率分布は、ｎ個作成される。

【0060】

＜同時確率分布：ステップＳ１０８＞
ｎ個のアミノ酸から任意の２つのアミノ酸を選び、一方の確率分布をｐ（Ａｘ）とし、他方の確率分布をｐ（Ａｙ）とする。計算時刻ｔの時のそれぞれの主要距離ＬをＬｘ（ｔ）、Ｌｙ（ｔ）とし、同一時刻同士のＬｘ（ｔ）とＬｙ（ｔ）を、ｘ軸ｙ軸で表される平面にプロットすることで、散乱分布を得ることができる。

【0061】

これをステップＳ１０６の場合同様に、確率分布ｐ（Ａｘ，Ａｙ）に変換する。これは、ｐ（Ａｘ）とｐ（Ａｙ）の同時確率分布を求めたものと考えることができる。また、これはＬｘ（ｔ）とＬｙ（ｔ）から求められた確率分布ｐ（Ａｘ）とｐ（Ａｙ）を周辺確率分布とする同時確率分布を求めたともいえる。同時確率分布ｐ（Ａｘ，Ａｙ）は、２つのアミノ酸から作成される。つまり、ｎ個のアミノ酸からはｎＣ２個の同時確率分布ｐ（Ａｘ，Ａｙ）が算出されることになる。

【0062】

＜相互情報量：ステップＳ１１０＞
相互情報量を図２（ｄ）のように求める。

【0063】

＜相互情報量マップ：ステップＳ１１２＞
相互情報量マップを図２（ｅ）のように、アミノ酸同士の相互情報量を表すように出力し、処理は終了する（ステップＳ１１４）。

【0064】

以上のようにして、原子間距離だけを用いて、平衡状態における各アミノ酸同士の動的相関量を求め、表示することができる。なお、本発明においては、上記の各ステップを実行するプログラムが含まれる。また、上記各ステップを実行する機能を実現する計算装置も本発明に含まれる。

【実施例0065】

ＰＤＺドメインタンパク質のアミノ酸間相互情報量を計算した結果を示す。ＰＤＺドメインは、様々な足場タンパク質において共通してみられる、８０から９０アミノ酸からなるタンパク質ドメインである。ＰＤＺドメインを有するタンパク質はＰＤＺドメインタンパク質と呼ばれている。

【0066】

ＰＤＺドメインは、リガンドタンパク質のＣ末端に結合する。この結合によって、シグナル伝達や細胞内輸送に関わる大きなタンパク質複合体が形成される。したがって、ＰＤＺドメインの動的な相関を知ることは重要である。

【0067】

また、ＰＤＺには、リガンドが結合した構造と、結合していない構造の２つが実験的に求められているので、それらの構造をシミュレーションに用いた。各々の構造はＸ線結晶構造解析で求められており、水素原子が欠損しているので水素原子を含む適正な分子となるようにこれを公知の方法で補完した。水素原子が補完された構造をモデル水分子で作られた直方体の中心に置き、タンパク質構造と重なり合う水分子を取り除いた。またこの直方体内の電荷がゼロになりかつ１５０ｍＭの濃度になるようにＮａ^＋イオンとＣｌ^－イオンを配置し重なり合う水分子を排除しシミュレーションに用いる系を構築した。この分子群で構成された直方体を周期的境界条件に置き分子間力を計算し、全系の原子に対するニュートン方程式を差分的に解き全ての分子内の原子を時間発展させ、所定の時間長に達するまで計算した。またステップＳ２１０で求める主成分軸は、共分散行列の固有値の占有率が８０％以上になるように求めた。

【0068】

図１０にＰＤＺドメインの立体構造をリボンモデルで示す。図１０（ａ）は、ＰＤＺドメインのみのリボンモデルである。図１０（ｂ）は、ＰＤＺドメインにリガンドが結合した図である。

【0069】

図１１（ａ）はリガンドペプチド分子が結合していない場合（図１０（ａ）の場合）の、アミノ酸側鎖間の相互情報量を示し、図１１（ｂ）はリガンドペプチド分子が結合した場合（図１０（ｂ）の場合）のアミノ酸側鎖間の相互情報量を示している。縦軸横軸はそれぞれＰＤＺドメイン内のアミノ酸番号である。また、グラフ右端に示すように、相互情報量は濃淡で表される。なお、相互情報量は、同じものの相互情報量は意味がないので、グラフの左下から右上に向かう対角線上は強制的に０（ゼロ）としている。

【0070】

図１１を見るとリガンドペプチドが結合していなかった際に、互いに動的相関量が多かったＰＤＺａや、ＰＤＺｂといった箇所にリガンドペプチドが結合すると、アミノ酸間の相関が有意に抑制されている事が観測できる。

【0071】

この図から明らかなように、本発明によってＰＤＺドメインにリガンドが結合した際のアミノ酸側鎖間の相関変化が大変明確に示されている。すなわち、本発明によってタンパク質等の生体高分子に関する分子内相関量が精度よく評価できることが解る。

【産業上の利用可能性】

【0072】

本発明は、シミュレーションを用いて分子内及び分子間の動的相関量解析を行う際に好適に利用することができる。

【符号の説明】

【0073】

１ コンピュータ
１０ 中央処理装置
１２ メモリ
１４ 外部メモリ
１６ 入力デバイス
１８ディスプレイ

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】