QM/MM法に基づいた溶液系のプロトン移動シミュレーション

Abstract

水素イオンは生体分子の理解に重要だが，分子シミュレーションにおいて最も取り扱いが難しい対象で，従来小さな系でのダイナミクスしか再現できなかった．今回，ハイブリッド分子モデルQM/MM法を拡張することでバルクでのプロトン輸送を再現した．これによりタンパク質など巨大系のプロトン輸送も取り扱い可能になる．

1.  はじめに

水素イオン（プロトン）は我々にとって馴染みの深い存在であり，化学・工学的そして生物学的にも重要である．この水素イオンの移動モデルは200年前にvon Grotthussによって提唱された．一方で水溶液中の水素イオンの安定構造については実験的な同定が難しく，Eigenカチオン（H₃O⁺）かZundel（H₅O₂⁺）かで長らく議論が続いた（図1）．この問題の解決には，計算機の発展を背景に1990年代に行われたTuckermannらの分子シミュレーションによる研究が大きく寄与した¹⁾．一般的に分子シミュレーションにおいて水素イオンは最も取り扱いが難しい存在である．それはGrotthuss機構で説明されるようにプロトンが周囲の水分子との共有結合の生成・消滅によって移動していくことにある．このことは電子構造を取り扱うことが可能な量子化学的（Quantum Mechanical: QM）モデルの必要性を意味する．このQMモデルを用いた分子動力学シミュレーションはab initio Molecular Dynamics（AIMD）と呼ばれる．AIMDは計算コスト（計算時間）が非常に大きく，モデルにもよるが対象系の大きさ（粒子数N）に応じてN ³からexp(N)のオーダーで増加する．したがって初期の水素イオンのAIMDシミュレーションも水32分子程度の系であった（これでは本誌読者の興味があろう生体分子系からは程遠い）．当時から計算機は大きく発展を遂げているが，AIMDの実行可能な対象系の大きさが劇的に変化しているわけではない．むしろ計算手法の工夫により大幅なコスト削減を図るアプローチが行われている．分割統治法に基づくアプローチがその最たる例であろう．これは系全体をいくつかのサブシステムに分割して個別にQMモデルを適用し最終的に統合する方法である，近年はバクテリオロドプシンなど生体分子への適用事例も報告されている²⁾．この方法では系のサイズに対するコストの増加は最大でN ^1.2程度に抑制されるが，それでも京／富岳コンピュータなどの数万コア以上の並列を前提としており，ラボレベルの計算機で気軽に実行することはできない．

図1

Eigenカチオン（H₃O⁺）（左上）とZundelカチオン（H₅O₂⁺）（左下）．Grotthuss機構の模式図（右図）．水素原子を白，酸素原子を赤で表示．

2.  QM/MM分子モデルとその問題点

量子と古典（Molecular Mechanics: MM）のハイブリッド分子モデルであるQM/MM法を用いたアプローチの開発も進んでいる．QM/MMモデルはタンパク質の活性部位など電子構造が重要な局所的な部位をQMモデルで計算し，それ以外の周辺環境をMMモデルで粗く計算するモデルである．この方法では系全体のサイズが計算の律速にはならないので，小／中規模の計算機でタンパク質など巨大分子の反応解析をする際に広く用いられてきた．このQM/MM法をプロトン移動の計算に利用すると考えるのは自然な発想であるが，プロトン移動ならではの2つの問題がそれを阻んできた．

1つ目の問題は溶媒の取り扱いである．水素イオンは溶媒の水分子と余剰水素の受け渡しをする．つまり溶質だけでなく溶媒もQMモデルで取り扱う必要がある．しかし溶質を取り囲むいくつかの水分子をQMモデルで定義しても，時間発展計算を行えばQM水分子は拡散し，代わりにMMの水分子が水素イオンを取り囲む．結果として溶質と溶媒の相互作用は古典的なものに置き換わるので水素イオンは取り扱えない．この問題の原因は，従来のQM/MM法ではMD計算の初めに分子毎にモデルを決めたら，計算中には切り替えることができないことにある．金属イオンの溶媒和など水の電子構造が重要な系でも同様の問題が発生する．

水素イオンに対する2つ目の課題は，溶質そのものの取り扱いである．たとえ後述するように溶質との距離により溶媒のモデルを切り替えることで，溶質をQM溶媒が取り囲むようにできても，水素イオンにおいては溶質そのものが時事刻々と変化する．したがって溶質とみなす対象も計算の最中で切り替わるような枠組みが必要となる．そのためにはまず余剰プロトンの座標を数値化する必要がある．この数値指標はExcess Proton Indicator（EPI）と呼ばれる．過去にもEPIの定式化が試みられたが，以下のような問題が生じていた（1）不連続性：ポテンシャル局面が不連続になる，（2）不安定性：プロトン移動の際にEPIが高速に変化し計算が不安定になる，（3）サイズ依存性：系が大きくなると水素イオンの座標を正しく示せない（4）高計算コスト．過去に提案されたEPIは（1, 2）を回避しているが（3, 4）の問題が顕在化するもの³⁾あるいはその逆の例はあったが⁴⁾，4つのすべての問題に対処できたEPIは存在しなかった．

3.  Full adaptive QM/MM法

前者の溶媒問題への対策としては，（solvent）adaptive QM/MM法が提案されてきた．これは溶質との距離に応じてMD計算の最中に溶媒の分子モデルが量子と古典の間で滑らかにスイッチするというものである．

コンセプトはシンプルだが，実際に空間的／時間的の連続性の問題が生じるため実現は容易ではなく，様々なadaptive QM/MM法が提案されてきた^5)-7)．しかし各手法それぞれ問題を抱えており解決には至らなかった⁸⁾．そこで我々は2014年にsize-consistent multipartitioning（SCMP）法を提案した⁹⁾．この方法は時間的な不連続性を除去し，現在までのところ最も空間的不連続性を抑制することができる方法であり，より正確に溶媒の量子化学効果をダイナミクスに取り込むことが可能である．

これに加え，今回我々は前述の4つの問題を回避できるようなEPIを提案した．まず溶質の問題（1）と（2）を回避するためにEPIに求められる条件は「系に含まれる水素イオンと水分子中のすべての酸素と水素原子の座標を含む連続関数」である．これを満たしつつ問題（3）と（4）を回避するために「水分子と水素イオンの局所的な溶媒和構造の違いを定量化し抽出できるような関数」となるようにEPIに工夫を加えた．さらにこのEPIをSCMP法に組み込むことで，full adaptive QM/MM法（溶質と溶媒の両方をadaptiveに取り扱える方法という意）として提唱した．そして実際にバルクの水溶液系でテストを行いプロトン移動のシミュレーションに成功した¹⁰⁾．

これにより様々な水素イオンのダイナミクスや物理化学量の定量的な解析が現実的な計算コストで実行が可能になる．実際，図2に示すように系のサイズが2倍，4倍と変化しても計算にかかる時間はほとんど変化しないことが確認できる．これは今回の手法がQM/MMに基づいていることによる利点である．さらに今回はバルク水中の水素イオンとして，水2047分子と水素イオンの系を用いて，水素イオンのダイナミクスを解析した．この際100 psのシミュレーションを独立して何回か実行して統計量を見積もった．

図2

系のサイズとMD 1 stepあたりの計算時間．

その結果の一部を以下に紹介する．図3に示すように真空中ではZundelカチオンが安定である．一方，溶液中ではEigenカチオンが逆転して安定になる．つまり溶液中では，主にEigenカチオンにあるプロトンがZundelカチオン構造を経由して隣の水分子へ輸送されていくことが定量的に示されている．さらに安定した長時間の計算は3,000回以上プロトン輸送イベントを平均化して時系列で表すことも可能にした．図4はプロトン輸送に伴ったバリア越えに際して，酸素間距離と酸素―水素間距離の振動がカップリングしている様子も鮮明に示している．この他にも溶媒和構造や拡散係数なども定量的な解析などが可能になったことが確認された．

図3

真空中でのH₃O⁺とH₂Oの2量体のポテンシャルエネルギー（上）と溶液中の自由エネルギー局面（下）．横軸は水素イオンと最近接する水分子の酸素間距離，縦軸は輸送される水素原子の酸素2原子の中心からのズレを表す．

図4

プロトン輸送の際の酸素原子間距離（黒線）と輸送されるプロトンの酸素2原子中心からのズレ（赤線）．対象の水素が酸素の中心に来た時刻をt = 0として表示．

4.  おわりに

今回提案したfull adaptive QM/MM法は，計算コストを抑制できるとはいえ，今後も利用できるのはたかだか数nsのMD計算となりそうである．したがって，生体分子系の研究に活用していくためには，より効率的なサンプリングが必要不可欠である．幸いにも今回の方法を用いれば水素イオンに外から人為的な力を加えることも可能になる．これは今回提案した方法が，外部力を加えてサンプリング効率を高めるアンブレラサンプリングやメタダイナミクス法などの計算技法との相性が良いことを意味しており，更なる展開の可能性を秘めている．

水素イオンは生体分子の構造や機能を議論する上で欠かせない存在である．しかしながら肝心の水素イオンを扱えないことは，分子シミュレーションにとって致命的な問題であり，筆者は学生の頃からその解決を目標にしてきた．同手法にはQM/MM法ならではのアーティファクトなど解決すべき課題は残っているが⁸⁾，それらの改善と同時に今後は様々な生体分子の機能の解析にも役立てていきたい．

最後に一連の研究に協力／サポートしてくださった共同研究者と大学関係者，およびJSTさきがけおよび日本学術振興会に感謝申し上げたい．

文献

• 1)   Tuckerman,  M. et al. (1995) J. Phys. Chem. 99, 5749-5752. DOI: 10.1021/j100016a003.
• 2)   Ono,  J. et al. (2020) J. Phys. Chem. B 124, 8524-8539. DOI: 10.1021/acs.jpcb.0c05507.
• 3)   Chakrabarti,  N. et al. (2014) Structure 12, 65-74. DOI: 10.1016/j.str.2003.11.017.
• 4)   Pezeshki,  S.,  Lin,  H. (2015) J. Chem. Theory Comput. 11, 2398-2411. DOI: 10.1021/ct501019y.
• 5)   Heyden,  A. et al. (2007) J. Phys. Chem. B 111, 2231-2241. DOI: 10.1021/jp0673617.
• 6)   Bulo,  R. et al. (2009) J. Chem. Theory, Comput. 5, 2212-2221. DOI: 10.1021/ct900148e.
• 7)   Takenaka,  N. et al. (2012) Chem. Phys. Lett. 524, 56-61. DOI: 10.1016/j.cplett.2011.12.053.
• 8)   Watanabe,  H. C.,  Cui,  Q. (2019) J. Chem. Theory. Comput. 15, 3917-3928. DOI: 10.1021/acs.jctc.9b00180.
• 9)   Watanabe,  H. C. et al. (2014) J. Chem. Theory. Comput. 10, 4242-4252. DOI: 10.1021/ct5005593.
• 10)   Watanabe,  H. C. et al. (2020) Phys. Chem. Chem. Phys. 23, 8344-8360. DOI: 10.1039/D1CP00116G.

Biographies

渡邉宙志（わたなべ　ひろし）

慶應義塾大学量子コンピューティングセンター特任講師