非調和性と環境を考慮した生体分子の振動計算

Abstract

振動分光法は分子の構造変化を鋭敏にプローブできる有力な解析法だが，スペクトルの帰属と解釈が困難であった．近年，分子の非調和性・環境・構造ゆらぎを考慮することで，生体分子の振動スペクトルが計算可能になった．その基盤技術である分子動力学計算，局所振動モードに基づく振動計算，QM/MM計算を紹介する．

1.  はじめに

振動分光法は，分子に光を照射したとき，分子の振動運動に共鳴する周波数において，光の吸収や散乱が起こることを利用し，分子の構造や機能を調べる手法である．特に，振動スペクトルは水素結合などの分子間相互作用を鋭敏にプローブできるのが特徴で，様々な生体分子の解析に用いられている．例えば，カルボニル基が水素結合しているときとそうでないとき，C-O結合距離はほとんど変わらない（<0.02 Å）が，C-O伸縮振動は50 cm^–1以上の有意な変化が観測される．このような特徴を活かし，アミド基に由来する振動バンドを解析することで，蛋白質の構造変化と機能が研究されている^1),2)．しかし，一般に生体分子の振動スペクトルは複雑で，その解釈は容易でない．そのため，理論計算による支援が重要である．

分子振動理論の出発点は調和近似である．調和近似は化学結合をバネに見立てる近似であり，ばね定数k（ポテンシャルの2次微分）から振動スペクトルを簡便に計算できるのが利点である．しかし，調和近似の問題点は，本物のポテンシャルが調和ポテンシャルと異なることである．この両者の差を「非調和性」と呼ぶ（図1a）．振動スペクトルを精度良く計算するには非調和性の考慮が必要不可欠である³⁾．

図1

（a）調和近似と非調和ポテンシャルの違い．ω，x_eは調和振動数と非調和定数．（b）QM/MM法の模式図．

さらに，生体分子を対象とする計算では，着目する分子振動を与える化学種が置かれた周囲環境も考慮する必要がある．振動計算はポテンシャルの精度に大きく影響される．そのため，ポテンシャルは高精度な量子化学計算で求める必要がある．しかし，量子化学計算は原子数とともに計算負荷が莫大になる．そこで，対象とする化学種のみを量子化学（QM）計算で扱い，周囲の蛋白質や溶媒などの環境を比較的計算負荷が低い分子力場（MM）で扱うQM/MM法（図1b）を利用する．QM/MM計算と非調和振動計算を組み合わせることで，生体環境にある化学種の振動スペクトルを高い精度で計算可能となった^4)-6)．

我々は分子動力学（MD）計算プログラムGENESIS^7),8)にQM/MM計算を実装し，これを振動計算プログラムSINDO⁹⁾と組み合わせた．GENESISとSINDOはいずれも我々が独自に開発したプログラムである．本稿では，これらのプログラムの実装を解説し，主要な成果として，バクテリオロドプシンの内部水に対する振動計算の結果を紹介する．

2.  GENESISにおけるQM/MM計算の実装

GENESISは理研を中心に開発したMD計算プログラムで，フリーソフトとして以下のURLから公開している（https://www.r-ccs.riken.jp/labs/cbrt/）^7),8)．GENESISの特徴は，高い並列性能を持つことと階層的手法（粗視化・全原子・QM/MM）が実装されていることである．全原子MD計算は，CHARMM力場とAMBER力場に対応し，「富岳」などの超並列スパコンを用いた高速計算が可能である．粗視化MDには，Goモデルやその拡張であるAICG2+モデルなどが利用可能である．また，レプリカ交換分子動力学法を始めとする，様々な拡張アンサンブル法と自由エネルギー計算法が実装されている．MD計算の結果（トラジェクトリー）を解析するプログラムも付随している．上記URLからアクセスできるチュートリアルには計算法の解説と実例が記述されている．

GENESISのQM/MM計算には2つの特徴がある．1つは高いポータビリティーである⁴⁾．QM/MM計算では，外部のQMプログラムを呼び出し，現在の座標情報を与えて，QMエネルギーとその座標微分など，必要な情報を受け取る．GENESISでは，外部QMプログラムの入力部分の相違を吸収できるインターフェースを個別に用意することで，様々なQMプログラムと容易に接続できる．表1に利用可能なQMプログラムを挙げる．GENESISのインプットには，利用するQMプログラム名，QMプログラムを実行するスクリプトの名前，QM計算のテンプレートファイルの名前を指定する．QM計算のオプションはテンプレートファイルに指定する．GENESISのインプットオプションを最小限にすることで，分かりやすく間違いを防ぐ仕様となっている．この他，SCF計算を前ステップの分子軌道からリスタートする，SCF計算が収束しなかったときに（QMオプションを変更して）再計算するなど，QM/MM-MD計算を実行する上で有用な様々な機能が実装されている．

表1

GENESISで利用可能なQMプログラム．

名前	備考
Gaussian	有償．最も広く使われている．
Q-Chem	有償．
DFTB+	フリー．DFTB計算が利用可能．
ORCA	アカデミアフリー．
TeraChem	有償．GPUによる高速計算可能．
QSimulate-QM	有償．並列による高速計算可能．

もう1つの特徴は，QSimulate-QM を利用した高速計算である．QSimulate-QMは並列計算を念頭にスクラッチから作り直されたプログラムである．GENESISとQSimulate-QMを並列計算機で実行することで，極めて高速な計算が可能である．密度汎関数（DFT）計算では10 ps/day，density functional based tight-binding（DFTB）計算では1 ns/dayのパフォーマンスが得られる（文献5の図1参照）．GENESISとQSimulate-QMは動的ライブラリーで接続されている．ファイルを介さないため，ディスクI/Oによる遅延を受けない．この仕様はDFTBによる高速計算において重要である．

また，反応経路探索法の1つであるString法を実装した⁵⁾．これを用いることで，酵素反応などの化学反応に応用可能である．

3.  グループ局所化座標に基づく振動計算

分子振動計算の目的は振動Schrödinger方程式を解くことである．その解法として，平均場近似であるvibrational self-consistent field（VSCF）法を出発点として，様々な拡張が提案されている．我々は擬縮退摂動論を用いたvibrational quasi-degenerate perturbation theory（VQDPT）法を開発し，SINDOプログラムに実装した．理論³⁾と実装⁹⁾の詳細は文献を参照されたい．

VQDPT法は分子サイズが100原子程度まで扱うことができるが，より大きな生体分子へ応用するには，さらなる方法論の拡張が必要である．ここで重要なのは，生体分子では，多くの場合，全ての振動モードを研究対象とするのではなく，生体機能に関連する部分構造に着目し，それに由来する振動モードの振動解析がなされることである．例えば，次節に示すバクテリオロドプシンでは，光受容プロトンポンプに重要な役割をする内部水の振動バンドが研究されてきた¹⁰⁾．このような背景から，我々はグループ局所化座標に基づく振動計算法を提案した⁶⁾．この方法では，振動座標を局所化させる分子，置換基，残基などを選ぶ．これをグループと呼ぶ．グループに局所化した振動座標に基づき振動ハミルトニアンを構築する．このとき，グループ内は非調和性を取り込むが，グループ間は調和カップリングだけに留める．すなわち，

  

$${\widehat{H}}_{vib}=\sum _g{\widehat{H}}_{vib}^{\left(g\right)}+\sum _{g,g\text{'}}{c}_{g_i{g}_j^{\text{'}}}{Q}_i^{\left(g\right)}{Q}_j^{\left(g\text{'}\right)},$$（1）

  

$${\widehat{H}}_{vib}^{\left(g\right)}=-\frac{1}{2}\sum _i\frac{{\partial }^2}{\partial {Q_i^{\left(g\right)}}^2}+V\left(Q_1^{\left(g\right)},\cdots ,Q_f^{\left(g\right)}\right).$$（2）

ここで $g$ はグループに対するインデックスで， $Q_i^{\left(g\right)}$ はグループ $g$ のi番目の基準座標である．

この方法では，グループ間の非調和性を無視するため，計算負荷を大幅に減らすことができるのがメリットである．一方，グループの定義はユーザの手に委ねられている．水分子や残基の側鎖など，直感的に分割できるときは良いが，そうでないと使えない．それでも，生体機能に関わる分子群が特定されている場合，このような近似法は極めて有用である．例として，膜蛋白質の内部水に着目した計算例を次節に紹介する．

4.  計算事例：バクテリオロドプシンの内部水

バクテリオロドプシンは光を受容し，プロトンを細胞内から外へ濃度勾配に逆らって輸送するプロトンポンプである．その光受容体はレチナール補因子である．レチナールは光受容によりtrans-cis光異性化反応を起こし，これがトリガーとなって，蛋白質の構造変化とプロトン移動反応が連続的に引き起こされる．レチナールの周辺には，プロトン化シッフ塩（PSB），極性残基と内部水により構成される水素結合ネットワークが発達している（図2a）．

図2

（a）バクテリオロドプシン全体と内部水周辺の拡大図（b）水素結合ネットワークの模式図．青数字はO-D，N-D伸縮振動モードの番号．（c）下から順に，BRに対する調和近似，非調和VQDPT2計算，構造平均を考慮した計算スペクトルと実験で観測された赤外差スペクトル．実験スペクトルでBRは負側に相当する．青い部分は強い水素結合に起因するバンド．Adapted with permission from [6]. Copyright 2021 American Chemical Society.

以前，柴田と神取¹⁰⁾は，重水においてresting state（BR）とレチナールの異性化反応直後（K state）の赤外差スペクトルを取得し，2050-2350 cm^–1に強い水素結合に由来するブロードな振動バンドを観測した（図2c）．詳細な解析から，この振動バンドは3つの内部水のO-D伸縮振動とPSBのN-D伸縮振動に帰属された（図2b, c）．このような強い水素結合は，プロトンポンプ機能を持つロドプシンに共通の特徴で，強い興味が持たれている．この振動バンドに対する従来の計算は調和近似に留まっており，非調和性と環境の影響を考慮した計算はなされていなかった．

QM/MM法と非調和振動計算を用いて，BR内部水の振動スペクトルを計算した．QM領域は，水素結合ネットワークを構成する分子（内部水，レチナール，R82，D85，D212）だけでは不十分で，それらに水素結合している残基（Y57, T89, Y185）と近傍の水分子（Wat403, 407）も含めた．QM原子数は137個で，計算レベルはB3LYP-D3/aug-cc-pVDZを用いた．また，非調和計算では，各水分子に局所化した座標を用いて，式（1）のように，水分子内の非調和性を考慮し，分子間は調和結合に留める近似を用いた．

図2cの下2つに，結晶構造から求めた調和および非調和計算のスペクトルを示す．両者の比較から，非調和性は極めて大きな影響があることが分かる．これは水素結合で一般的に見られる傾向である．一方，非調和計算において2000 cm^–1付近に得られる大きいピークが実験には存在しない．このピークは，モード47（Wat402の片方のO-D伸縮）に起因するが，実験的には2180 cm^–1のピークに帰属されており，大きく異なっている．

次に，BRを脂質二重膜に埋め込み，古典MD計算を実行した．得られたトラジェクトリーの中から8つのスナップショット構造を選び，同様にQM/MM非調和振動計算により振動スペクトルを計算した．それらを平均した結果を図2cの上から2番目に示す．スナップショットの選び方には注意を要するが（詳細は文献6），構造平均を取ることで実験を良く再現した．興味深いことに，振動バンドはそれぞれ1つのO-D/N-D伸縮振動に対応するわけではなく，Wat402（モード47），Wat401（モード38），PSB（モード18）の重ね合わせであった．これは，PSB-Wat402-D85-Wat401をつなぐ水素結合が特に強く，それが動的に揺らぐことを反映している．

5.  今後の展望

局所化座標に基づく非調和振動計算とQM/MM計算を組み合わせることで，生体分子の振動解析が可能になった．計算と実験スペクトルを比較することで，計算方法を検証し，構造・機能に対する理解が得られた．今後は，さらにバクテリオロドプシンの光反応中間体や他のロドプシンへ展開する．また，アミドバンドの新しい計算法に取り組み，蛋白質の構造変化を振動スペクトルから検出する方法の開発に挑戦する．

文献

• 1)   Ganim,  Z. et al. (2008) Acc. Chem. Res. 41, 432-441. DOI: 10.1021/ar700188n.
• 2)   Feng,  C. -J. et al. (2021) J. Phys. Chem. B 125, 4620-4633. DOI: 10.1021/acs.jpcb.1c00399.
• 3)   Yagi,  K. (2016) Mol. Sci. 10, A0085. DOI: 10.3175/molsci.10.A0085.
• 4)   Yagi,  K. et al. (2019) J. Chem. Theory Comput. 15, 1924-1938. DOI: 10.1021/acs.jctc.8b01193.
• 5)   Yagi,  K. et al. (2021) J. Phys. Chem. B 125, 4701-4713. DOI: 10.1021/acs.jpcb.1c01862.
• 6)   Yagi,  K.,  Sugita,  Y. (2021) J. Chem. Theory Comput. 17, 5007-5020. DOI: 10.1021/acs.jctc.1c00060.
• 7)   Jung,  J. et al. (2015) Wiley Interdiscip. Rev. Comput. Mol. Sci. 5, 310-323. DOI: 10.1002/wcms.1220.
• 8)   Kobayashi,  C. et al. (2017) J. Comput. Chem. 38, 2193-2206. DOI: 10.1002/jcc.24874.
• 9)   Yagi,  K.（2022）アンサンブル 24, 122-129. DOI: 10.11436/mssj.24.122.
• 10)   Shibata,  M.,  Kandori,  H. (2005) Biochemistry 44, 7406-7413. DOI: 10.1021/bi050122+.

Biographies

八木　清（やぎ　きよし）

理化学研究所開拓研究本部専任研究員

杉田有治（すぎた　ゆうじ）

理化学研究所開拓研究本部主任研究員，同計算科学研究センターチームリーダー，同生命機能科学研究センターチームリーダー