Nuclear Quantum Effect on the Intermolecular Hydrogen Bond of Acetic Acid − phosphorous Acid Anion Cluster: an ab initio Path Integral Molecular Dynamics Study

Abstract

High Performance Computerで効率の高いシミュレーションを実行できる，分子科学計算ソフトウェアNTChemと経路積分分子動力学(PIMD)法を統合した階層的並列プログラムプラットフォームを用いることにより，低障壁水素結合(LBHB)の存在が示唆されているperiplasmic phosphate binding protein (PPBP)のモデル分子である酢酸－リン酸アニオンクラスターの分子間水素結合における原子核の量子揺らぎの効果を明らかにした．モデル分子における一つの水素結合を形成する配置を解析したところ，原子核を古典的な質点で取り扱った古典シミュレーションでは，プロトンは水素結合の両端に局在化し，通常の水素結合と同様の振る舞いを示した．一方，原子核の量子揺らぎの効果を取り込んだ量子シミュレーションでは，プロトンが水素結合の中央に存在する確率が高くなった．酢酸－リン酸アニオンクラスターの分子間水素結合は，LBHBと同じような傾向を示した．

1 研究目的

プロトン移動反応の障壁が低くなる，いわゆる低障壁水素結合(Low Barrier Hydrogen Bond (LBHB))は，近年多くの注目を集めており，実験と理論の両面から研究が進められている [1,2,3,4,5]．LBHBは，通常の水素結合よりも短く，結合が強くなることが知られている．また，通常の水素結合ではプロトンが水素結合の両端の重原子のいずれかに偏るのに対して，LBHBでは，プロトンは両端の重原子の中央に位置する．このLBHBが生体分子における触媒において重要な役割を果たすことが示唆されている [6,7,8,9,10,11,12]．最近では，山口らはX線結晶解析や中性子散乱実験にてLBHBの存在を見出した [9]．しかし，量子化学計算を用いた理論解析では，その存在を否定する結果が得られており [13]，さらなる解析が望まれている．

イオンを含む分子クラスターにおいては，通常の量子化学計算では安定な構造として見いだすことができないLBHBを正確に取り扱うためには，核の量子揺らぎと温度揺らぎを考慮することが必要であることが明らかにされつつある [14,15,16,17,18]．その双方を正確に取り扱うことが可能なab initio経路積分分子動力学(PIMD)法 [19,20,21,22]を用いた分子シミュレーションはLBHBを解析するための最良のツールの一つとなっている．

通常の電子状態計算では原子核を一つの微小な質点とする古典的なアプローチを用いるのに対して，PIMD法では，ばねで結合した多数の古典的な質点(ビーズ)の集合で原子核を表現する．分子シミュレーションで発生したそれぞれの古典的なビーズの経路のアンサンブル平均が原子核の量子揺らぎを表現している．PIMD法では，それぞれのビーズにおいてab initio電子状態計算の実行が必要となり，通常の電子状態およびそれに基づく分子シミュレーションよりも多くの計算コストを要する．特に，注目されている大規模生体分子におけるLBHBを取り扱うためには，High Performance Computer (HPC)の利用が欠かせない．そこで，本研究ではHPCを効率よくPIMD法を用いたシミュレーションを実行できる計算基盤を確立する．具体的にはHPCを用いた効率のよい電子状態計算を実現しているNTChem分子科学計算ソフトウェア [23,24]とPIMD法のプログラムを連携させた階層的並列プログラムプラットフォームを構築することにより，HPCにおける効率のよいPIMDシミュレーションを実現する．

構築した新しいプラットフォームを用いてリン酸-酢酸アニオンクラスターのPIMDシミュレーションを実行する．その分子系をFigure 1に示す．この分子系は，periplasmic phosphate binding protein (PPBP)のモデル分子である．分子内にLBHBが存在するとされているPPBPはリン酸とヒ酸を識別し，毒性の高いヒ酸を生体内から取り出す機能を有することが報告されている [25]．また，タンパク質と酸の間に水素結合を形成するが，その水素結合がLBHBであり，そのLBHBがリン酸とヒ酸を識別することが示唆されている．そこで，このタンパク質のモデル分子とこれらの酸を結合させたモデル分子系の酢酸－リン酸アニオンクラスターのPIMDシミュレーションを実行し，酢酸－リン酸間の分子間水素結合における原子核の量子揺らぎの効果を検証した．

Figure 1.

 The acetic acid (left hand side) − phosphorous acid (right) anion cluster

2 計算方法

ab initio PIMD法は，温度効果に加え，on-the-flyで電子状態を量子力学的に，さらには核自身の量子性も含んだ，全自由度の量子力学的取扱いを可能とする手法である．on-the-fly PIMD法では，核の量子性を古典粒子の集まり(ビーズ)として表現し，古典粒子からなる分子系の計算にはab initio分子軌道法を用いる．この2のab initio分子軌道法については高並列分子軌道法プログラムNTChemを用いる．

我々が実装したon-the-fly PIMDプログラムについては，MPIおよびOpenMPを用いたhybrid並列を導入した．こちらについては既に並列プログラムを実装済みであったが，HPC環境で効率のよい並列計算が実現できるよう，通信やI/Oの高度化を行った．さらに，NTChemと我々のプログラムをシームレスに接続することを可能にするインターフェースを作成し，HPCに最適な階層的並列プログラムプラットフォームを構築した．

その構築したプラットフォームを用いて，PPBPのモデル分子である酢酸－リン酸アニオンクラスターのPIMDシミュレーションを実行した．酢酸分子とリン酸分子が水素結合を形成するモデル分子を用いて通常の原子核の量子揺らぎを考慮しない古典的なシミュレーション(以降，古典シミュレーション)と量子揺らぎを取り込む量子的なPIMDシミュレーション(以降，量子シミュレーション)を実行し，二つの酸の違いと水素結合に及ぼす核の量子揺らぎの効果について調べた．

また，モデル分子においては，完全な生体環境を実現できないため，様々な初期配置(古典シミュレーションは10個，量子シミュレーションは6個)を用いてそれぞれシミュレーションを実行した．シミュレーションは0.1 fs刻みで運動方程式を解き，分子系の温度を300 Kに設定した．熱浴としてNose-Hoover Chain [17] (長さ4)を用いた．量子シミュレーションは5,000ステップ，古典シミュレーションは10,000ステップの平衡シミュレーションを実行した上で，量子シミュレーションはそれぞれ25,000ステップ，古典シミュレーションはそれぞれ30,000ステップのプロダクションランを実行した．また，量子シミュレーションのビーズ数を16とした．

シミュレーションの初期配置については量子シミュレーションの初期配置をFigure 2に示す．分子間で形成する水素結合におけるプロトンの位置の異なる配置を初期配置として用いている．リン酸の酸素に近い位置にプロトンが存在する配置を3個(古典シミュレーションならば5個)，酢酸の酸素に近い位置にプロトンが存在するものを3個(古典シミュレーションならば5個)とそれぞれが50%になるようにした．

Figure 2.

 Six initial structures for quantum simulations of the acetic acid − phosphorous acid anion cluster. (a) Initial structures with the proton in the hydrogen bond close to the oxygen atom of acetic acid. (b) Initial structures with the proton in the hydrogen bond close to the oxygen atom of phosphorous acid.

電子状態計算および電子状態に基づく力は，NTChemソフトウェアを用いてωB97-XD [26]で計算した．PとAsの基底関数には，LANL2DZ [27,28,29]にdiffused関数を加えたECP [30]を用い，それ以外の原子については6–31++G**を用いた．

複数のシミュレーションを同時に実行し，得られたシミュレーションの解析を行った．その結果，本研究で取り扱った分子系では着目していた水素結合以外にも水素結合を形成した．水素結合を二つ形成した配置と三つ形成した配置のスナップショットをFigure 3に示す．なお，水素結合は，重原子間が3.0 Å以内のものとした．このようなケースでは，通常よりも多くの分子構造を用いて解析を行う必要があるが，複数の初期配置を用いてシミュレーションを行うことにより，様々な水素結合構造を効率よくサンプルすることができた．

Figure 3.

 A representative snapshot of a configuration obtained from simulation: (a) structure forming two hydrogen bonds and (b) structure forming three hydrogen bonds.

3 結果と考察

3.1 複数の分子間水素結合の形成

酢酸－リン酸アニオンクラスターのab initio PIMD分子シミュレーションを実行し，その結果の解析を行った．

生体分子PPBPにおいてはFigure 1で示すように一つの分子間水素結合を形成する．本研究のシミュレーションで得られた全ての配置において分子間水素結合を形成することを確認した．しかし，孤立したモデル分子を用いたため，一つの分子間水素結合だけではなく，二つ及び三つの分子間水素結合を形成する配置も見られた．得られた配置のうち，水素結合を一つ，二つ，そして三つ形成する配置の割合，及び，それらの配置における平均水素結合距離を全体の平均と併せてTable 1に示す．

Table 1.  The average hydrogen bond distance (Angstrom) of the configurations obtained from the molecular simulation of the acetic acid − phosphorous acid anion cluster. The percentage of the configurations with one, two, and three hydrogen bonds (HB) and their average hydrogen bond distance are listed as well.

	Quantum	Classical
All Configurations
ROO	2.628	2.657
RδH	0.492	0.440
Configurations with 1 HB	28%	42%
ROO	2.494	2.557
RδH	-0.250	-0.426
Configurations with 2HBs	61%	42%
ROO	2.628	2.641
RδH	0.577	0.603
Configurations with 3HBs	11%	16%
ROO	2.734	2.760
RδH	0.805	0.862

まず，分子間でいくつの水素結合を形成するか解析を行った．古典シミュレーションでは，水素結合を一つ，二つ，そして三つ形成する配置の割合がそれぞれ42%，42%，16％であった．それに対して，量子シミュレーションでは一つ，二つ，そして三つ形成する配置の割合がそれぞれ28%，61%，11％であった．

古典シミュレーションでは，水素結合を一つ形成する配置と二つ形成する配置がほぼ同数であったのに対して，量子シミュレーションでは，一つ水素結合を形成する配置の割合が少ないことがわかった．古典シミュレーションと量子シミュレーションにおいて，水素結合形成に優位な差が見られた．

3.2 水素結合における平均結合距離の解析

次に，古典シミュレーションと量子シミュレーションにおける分子間水素結合の平均距離について解析を行った．まずは，水素結合における重原子間距離について着目したところ，古典シミュレーションの結果よりも量子シミュレーションにおける距離のほうが水素結合の数にかかわらず，短いことがわかる．特に，一つの水素結合を形成した配置においては，大きな核の量子揺らぎの効果が見られた．水素結合の数が増えると差が少なくなった．水素結合を一つ形成するときに，大きな量子効果が見られた．

水素結合の数と重原子間距離の関係をより明確に捉えるために，分布についても解析を行った．重原子間距離の分布をFigure 4に示す．なお，水素結合の数と重原子間距離の関係の傾向は古典シミュレーションと量子シミュレーションにおいて差は少なかったので，量子シミュレーションの結果のみ示す．こちらのプロットは，それぞれの配置におけるすべての水素結合の重原子間距離についてプロットしたものである．水素結合を二つ形成する構造については二つの水素結合の重原子間距離，三つ形成する構造についても三つすべての水素結合の重原子間距離をプロットしている．そのため，平均水素結合距離と併せて水素結合1つ1つをより明確に捉えることができる．

Figure 4.

 The probability distribution of the hydrogen bond distance of R_OO for the quantum simulation. The dark blue line, the green line, the light blue line, and the yellow line, represent the distribution of all configurations, configurations with one hydrogen bond, two hydrogen bonds, and three hydrogen bonds, respectively.

濃青線のサンプル全体の分布を見ると，水素結合距離は2.25 Åから3.0 Å以上まで幅広く分布していることがわかる．緑線の水素結合を1個のみ形成するサンプルの分布を見ると特に2.25 Åから2.7 Åの領域に多く分布する．それに対して，2–3個の水素結合を持つサンプル(淡青線，黄線)は幅広く分布し，特に黄線の3個水素結合を形成するサンプルは，2.4 Å以上で多く分布する．水素結合の水素結合の数が増加するほど水素結合の重原子間距離が長くなっている様子がわかる．

また，緑線の水素結合を1個のみ形成するサンプルの分布は一つのピークを持つのに対して，3個の水素結合を持つサンプル(黄線)は二つのピークを持つことがわかる．水素結合を2個持つサンプル(淡青線)は3個の場合ほどははっきりとはしていないが，やはり2つのピークを持つ．それは，すべての水素結合が同じ結合距離となるわけではなく，短めの水素結合と長めの水素結合を持つことを示している．しかし，短距離側のピークを比較すると，水素結合が増えるに従い，ピークは長距離側へとシフトする．つまり，水素結合の水素結合の数が増加するほど水素結合の重原子間距離が長くなっている傾向は変わらない．

水素結合における核の量子揺らぎの効果についてより詳細に調べるために，水素結合におけるプロトンの位置の解析を行った．そのために，水素結合におけるプロトンの座標であるFigure 1のR_δHを定義した．Table 1が示すように水素結合を一つ形成する配置では，R_δHは負の値を示すのに対して，水素結合の数が増加すると正に変化することがわかった．すなわち，水素結合が一つの場合は，酢酸の方に偏り，水素結合の数が増えるとリン酸の方に偏る．次に，古典シミュレーションと量子シミュレーションを比較した．水素結合を一つ形成する場合は，量子シミュレーションは水素結合の中央付近にプロトンがあるのに対して，古典シミュレーションでは，酢酸の方に大きく偏る．水素結合を増やすと，古典シミュレーションの方が量子シミュレーションよりもリン酸側に偏る．つまり，水素結合の数に関わらず，量子揺らぎの効果を取り込んだ場合の方が，プロトンは水素結合の中央に近く位置することがわかった．水素結合の数，そして，核の量子揺らぎの効果を取り込むことにより，プロトンの位置が変化する様子を捉えることができた．

3.3 水素結合におけるプロトン移動

プロトンの位置やプロトン移動の様子をより詳細に捉えるために，水素結合におけるプロトン座標R_δHについてさらに解析を行った．プロトン座標R_δHの確率分布をFigure 5に示す．なお，こちらの分布は配置におけるすべての水素結合におけるプロトン座標をプロットしている．水素結合を二つ形成する構造については二つの水素結合のプロトン座標，三つ形成する構造についても三つすべての水素結合のプロトン座標をプロットしている．そのため，平均水素結合距離と併せて水素結合1つ1つをより明確に捉えることができる．まず，濃青線の全体のサンプルに着目すると，全てのシミュレーションにおいてプラス領域とマイナス領域にそれぞれピークを持つことがわかる．量子シミュレーションの結果と古典シミュレーションの結果を比較すると，古典シミュレーションでは，0近傍で分布がほぼなく，水素結合の両端の重原子へ偏るのに対して，量子シミュレーションでは0でも分布が多く見られた．LBHBにおいては，古典シミュレーションでは水素結合におけるプロトンが水素結合の両端の重原子へ偏り，量子シミュレーションでは水素結合の中央にプロトンの分布が増加する傾向にあり [14,15,16,17,18]，この水素結合はLBHBに似た傾向を示している．

Figure 5.

 The probability distribution of the hydrogen bond distance R_OO: (a) quantum simulation and (b) classical simulation. The dark blue line, the green line, the light blue line, and the yellow line, represent the distribution of all configurations, configurations with one hydrogen bond, two hydrogen bonds, and three hydrogen bonds, respectively.

次に，緑線の水素結合を一つ形成する配置における分布を見ると，古典シミュレーションも量子シミュレーションもマイナスの分布では1個の水素結合を形成するサンプルが大きく寄与していることがわかる．また，緑線のプロットについて量子シミュレーションと古典シミュレーションを比較すると，量子シミュレーションはプラス領域にも分布があるのに対して，古典シミュレーションでは分布が見られない．古典シミュレーションでは，水素結合が1個形成される時，酢酸の酸素に偏り，ほとんど動かないのに対して，量子シミュレーションでは，プロトン分布は中央に大きくシフトし，水素結合の両端にプロトンが広く分布することがわかる．

一方で，水素結合の数が増加すると，いずれのシミュレーションにおいても水素結合が増えるに従い，プラス領域の分布が多くなる．これは，リン酸の方の酸素にプロトンが偏ることを意味する．この場合においても，古典シミュレーションでは水素がほとんど動かないのに対して，量子シミュレーションでは水素結合の両端にプロトンが広く分布することがわかる．

本研究ではモデル分子を用いたことにより，分子間水素結合を一つだけ形成する配置だけでなく，水素結合を複数形成する配置が見られた．水素結合が複数存在する配置では量子効果はさほど見られなかったが，水素結合を一つ形成する配置については，原子核の量子揺らぎの効果を取り込んだ量子シミュレーションで，古典シミュレーションでは水素結合においてプロトンが水素結合両端の重原子に偏るのに対して，プロトンが水素結合の中央近傍に存在する確率が高くなり，LBHBと同じ傾向を示すことが明らかになった．

4 結論

HPC環境で効率の高いab initio経路積分分子動力学シミュレーションを実現する階層的並列プログラムプラットフォームを用いて，低障壁水素結合の存在が報告されている生体分子PPBPのモデル分子である酢酸－リン酸アニオンクラスターのシミュレーションを行い，分子間水素結合における原子核の量子揺らぎの効果を解析した．モデル分子のシミュレーションでは，実際のタンパク質では報告されていない複数の水素結合を持つ配置が見られた．水素結合を一つ形成する配置を中心に水素結合の構造を解析したところ，古典シミュレーションでは水素結合のプロトンが水素結合の両端に偏るのに対して，量子シミュレーションでは水素結合中央近傍のプロトンの分布確率が大幅に高くなることがわかった．水素結合が1個の時，PPBPの水素結合はLBHBに近い傾向を示すことが明らかになった．

現在，リン酸をヒ酸に変えたよりPPBPに近いモデル分子のab initio PIMD分子シミュレーションを実行し，その結果を解析している．リン酸とヒ酸の違いについて明らかにすることができれば，PPBPがリン酸とヒ酸を識別するメカニズムにせまることができる．開発したプラットフォームを用いることにより，実際の生体環境を考慮した大規模分子系での検証も視野に入る．今後，これまでに分かっていなかった識別メカニズムの解明に取り組む予定である．

参考文献

• [1]   G. A. Jeffrey, An Introduction to Hydrogen Bonding, Oxford University Press (1997).
• [2]   A. Kohen, H. H. Limbach, Isotope Effects in Chemistry and Biology, CRC Press (2005).
• [3]   T. Steiner, Angew. Chem. Int. Ed., 41, 48 (2002).
• [4]   W. H. Robertson, M. A. Johnson, Annu. Rev. Phys. Chem., 54, 173 (2003).
• [5]   M. Kabeláč, P. Hobza, Chemistry, 7, 2067 (2001).
• [6]   J. A. Gerlt, P. G. Gassman, Biochemistry, 32, 11943 (1993).
• [7]   W. W. Cleland, M. M. Kreevoy, Science, 264, 1887 (1994).
• [8]   P. A. Frey, S. A. Whitt, J. B. Tobin, Science, 264, 1927 (1994).
• [9]   S. Yamaguchi, H. Kamikubo, K. Kurihara, R. Kuroki, N. Niimura, N. Shimizu, Y. Yamazaki, M. Kataoka, Proc. Natl. Acad. Sci. USA, 106, 440 (2009).
• [10]   S. O. Shan, S. Loh, D. Herschlag, Science, 272, 97 (1996).
• [11]   D. B. Northrop, Acc. Chem. Res., 34, 790 (2001).
• [12]   T. Ishida, Biochemistry, 45, 5413 (2006).
• [13]   K. Saito, H. Ishikita, Proc. Natl. Acad. Sci. USA, 109, 167 (2012).
• [14]   M. Tachikawa, M. Shiga, J. Am. Chem. Soc., 127, 11908 (2005).
• [15]   K. Suzuki, M. Shiga, M. Tachikawa, J. Chem. Phys., 129, 144310 (2008).
• [16]   Y. Kawashima, M. Tachikawa, Chem. Phys. Lett., 571, 23 (2013).
• [17]   G. J. Martyna, M. L. Klein, M. Tuckerman, J. Chem. Phys., 97, 2635 (1992).
• [18]   M. E. Tuckerman, D. Marx, M. L. Klein, M. Parrinello, Science, 275, 817 (1997).
• [19]   D. Marx, M. Parrinello, Z. Phys. B Condens. Matter, 95, 143 (1994).
• [20]   D. Marx, M. Parrinello, J. Chem. Phys., 104, 4077 (1996).
• [21]   M. Shiga, M. Tachikawa, S. Miura, Chem. Phys. Lett., 332, 396 (2000).
• [22]   M. Shiga, M. Tachikawa, S. Miura, J. Chem. Phys., 115, 9149 (2001).
• [23]   T. Nakajima, M. Katouda, M. Kamiya, Y. Nakatsuka, Int. J. Quantum Chem., 115, 349 (2015).
• [24]   NTChem, 2013.http://labs.aics.riken.jp/nakajimat_top/ntchem_j.html
• [25]   M. Elias, A. Wellner, K. Goldin-Azulay, E. Chabriere, J. A. Vorholt, T. J. Erb, D. S. Tawfik, Nature, 491, 134 (2012).
• [26]   J.-D. Chai, M. Head-Gordon, J. Chem. Phys., 128, 084106 (2008).
• [27]   P. J. Hay, W. R. Wadt, J. Chem. Phys., 82, 270 (1985).
• [28]   W. R. Wadt, P. J. Hay, J. Chem. Phys., 82, 284 (1985).
• [29]   P. J. Hay, W. R. Wadt, J. Chem. Phys., 82, 299 (1985).
• [30]   C. E. Check, T. O. Faust, J. M. Bailey, B. J. Wright, T. M. Gilbert, L. S. Sunderlin, J. Phys. Chem. A, 105, 8111 (2001).