Development of Machine Learning Models withFragment Molecular Orbital Calculation Data

Abstract

フラグメント分子軌道(FMO)法はタンパク質全体を量子化学計算可能な稀有な手法である．そして，FMO法によって得られるデータもまた，現状ではタンパク質系の量子化学計算データとして唯一無二のものとなっている．汎用ソフトウェアでは生成が困難なタンパク質の量子化学計算データとそれを用いた様々な機械学習モデルの開発は，近年活性化が著しいAI創薬に大きなインパクトを与えることが期待される．本稿では，筆者らのグループで進めているFMOデータを用いた機械学習モデル(原子電荷予測モデル，相互作用予測モデル，機械学習力場)の開発状況を概説する．

Translated Abstract

Abstract: Fragment Molecular Orbital (FMO) is a unique method that allows quantum mechanical (QM) calculations of entire proteins. The data obtained by the FMO method are also currently the only QM calculation data for protein systems. The development of various machine learning models using the QM calculation data of proteins, which are difficult to generate with general-purpose software, is expected to have a significant impact on AI drug discovery, which has been remarkably active in recent years. This paper outlines the status of the development of machine learning models using FMO data, which is ongoing in the author's group.

1 はじめに

計算機ハードウェアの進化に伴い深層学習や強化学習が現実のものとなり，科学第4のパラダイムとしてデータ科学は広く認知されるに至った．2016年にGoogle DeepMind社の囲碁プログラムAlphaGoが当時のトッププロ棋士を破ってから早くも10年が経とうとしている．材料開発へのデータ科学の応用であるマテリアルズ・インフォマティクスも研究がますます活発化している．創薬分野においては，古くからケモインフォマティクスとしてデータ科学手法が活用されてきたが，近年ではAI創薬として新たなステージへと進化していると言っても過言ではない．これらデータ科学の活用は，材料探索・分子生成・実験条件最適化など多岐にわたるが，さらにシミュレーションを高度化・高速化するような取り組みも盛んに行われている．その最たるものは機械学習力場であろう．古典力場においては式(1)のようなポテンシャル関数によって原子間相互作用を表現 [1,2,3,4]する．ポテンシャル関数は基本的に結合長，結合角，二面角に関係する結合エネルギー項と，ファンデールワールス力及びクーロン力に起因する非結合エネルギー項より構成される．

  

$E_{total}=\sum _{bonds}{K}_r{\left(r-r_{eq}\right)}^2+\sum _{angles}{K}_{\theta }{\left(\theta -{\theta }_{eq}\right)}^2+\sum _{dihed}\frac{V_n}{2}\left[1+cos\left(n\varnothing -\gamma \right)\right]+\sum _{i<j}\left[\frac{A_{ij}}{R_{ij}^{12}}-\frac{B_{ij}}{R_{ij}^6}+\frac{q_i{q}_j}{\epsilon R_{ij}}\right]$(1)

この式(1)の右辺(及び定式化されていな高次の相互作用も含む)を全て機械学習モデルで代替するものが機械学習力場である．言い換えると，与えられた座標情報に基づいて，ポテンシャルエネルギーと個々の原子に働く力を予測する機械学習モデルとも言える．原理的には，学習データで加味されている相互作用は全て学習・予測可能である．機械学習力場には，ニューラルネットワーク(NN) [5]，グラフニューラルネットワーク(GNN) [6, 7]，ガウス過程回帰 [8]を用いるものなどが報告されているが，現在の主流はNNやGNNを用いたものである．機械学習力場の学習においては，第一原理計算・量子化学計算により生成されたデータセットが必要とされ，無機結晶材料系ではVASP(Vienna Ab initio Simulation Package) [9,10,11]やQuantum-ESPRESSO [12, 13]などのコードが用いられ，低分子系ではgaussian [14]によるデータ生成が一般的である．機械学習力場の学習プログラム・スクリプトもオープンソースとして公開されているもの [15, 16]が複数あり，機械学習力場の構築・適用は世界的な大きなブームとなっている．しかしながら，これら機械学習力場の生体高分子系への適用は極めて限定的 [17, 18]である．その理由として，生体高分子系の第一原理計算・量子化学計算による学習データ生成が極めて困難である点が挙げられる．Unkeらは，小タンパク質の切り出し構造に対して第一原理計算を行うことでデータを生成して対応 [18]した．一方，筆者らのグループはフラグメント分子軌道(FMO)法 [19]によりデータ生成の困難を克服している．

FMO法は1999年に北浦らによって考案された分子の電子状態の近似計算法 [19]である．その名のとおり計算対象をフラグメントに分割して計算を行うが，多体効果を効果的に取り込んでいることが大きな特徴であり，通常の量子化学計算と同等の結果を得ることができる．さらに，フラグメントに分割して計算するためフラグメント間相互作用(IFIE)を計算することができ，タンパク質‐リガンド複合体におけるタンパク質の個々のアミノ酸残基とリガンドの相互作用を，電子状態に基づいて定量的に解析することが可能である．さらに，原子電荷のような個々の原子特性を計算することも可能である．したがって，FMO法を使うことで①タンパク質単位②フラグメント単位③原子単位という3階層の量子化学計算データを同時に生成可能となる．筆者らのグループでは，この3階層のデータと機械学習力場のフレームワークを用いた機械学習モデルの開発を進めている．本稿ではそれらの開発状況について紹介する．

2 機械学習力場のフレームワーク

筆者らのグループでは，BehlerとParrinelloにより提案されたHDNNP(High-Dimensional Neural Network Potential) [5]をベースに開発を進めている．HDNNPでは元素毎にNNを構築する．系に含まれる各原子に対して，対応する元素用のNNを用いて原子エネルギーを算出し，それらの和(E_NN)によって系の全エネルギーを表現する．FMO法で得られた全エネルギー(E_FMO)とE_NNとの比較によって，各NNのパラメータを最適化する．各原子に働く力については，学習により得られたHDNNPを各原子座標に関して解析的に微分することで算出可能である．

元素ごとのNNに入力する記述子は，Behlerらにより提案されたACSFs(Atom Centered Symmetry Functions) [20]を主に用いている．ACSFsは，着目原子iの周辺環境を式(2)~(4)に示す関数により記述する．

  

$G_i^{\text{rad}}=\sum _{j\ne i}{e}^{-\eta {\left(r_{ij}-r_s\right)}^2}{f}_{c\left(r_{ij}\right)}$(2)

  

$G_i^{ang}=2^{1-\zeta }\sum _{j\ne i}\sum _{k\ne i,j}{\left(1+\lambda \cos \left({\alpha }_{ijk}\right)\right)}^{\zeta }{e}^{-\eta {\left(r_{ij}+r_{ik}+r_{jk}\right)}^2}{f}_{c\left(r_{ij}\right)}{f}_{c\left(r_{ik}\right)}{f}_{c\left(r_{jk}\right)}$(3)

  

$f_c\left(r_{ij}\right)=\{\begin{array}{c}tan{h}^3\left[1-\frac{r_{ij}}{r_c}\right]with{r}_{ij}\le r_c\\ 0with{r}_{ij}>r_c\end{array}$(4)

式(2)は着目原子iの周辺環境を動径分布として，式(3)は角度分布として記述する．また，これらの分布は，式(4)のカットオフ関数により着目原子iからカットオフ距離(r_c)以内に存在する原子についての寄与のみを考えるACSFsは分子の並進や回転，また原子の交換に対しても不変である．さらには，原子座標に関して解析的に微分可能であることから，機械学習力場に必要な各原子に働く力を計算することができるため，NNP構築には好適な記述子である．ACSFsの他にも，wACSFs(weighted-ACSFs) [21]やChebyshev多項式 [22]を用いるものも存在する．

これらHDNNPとACSFsを用いることで，FMO法によって得られる唯一無二の①タンパク質単位②フラグメント単位③原子単位の量子化学データの学習が可能である．具体的には，①のタンパク質単位のデータについては，機械学習力場のフレームワークをそのまま用いることで生体高分子用の機械学習力場の構築が可能と考えられる．②のフラグメント単位のデータについては，各元素のNN出力を原子IFIEとみなし，対象とするフラグメントを構成する原子の分だけ原子IFIEを足し合わせた上で，FMO法で得られた教師IFIEとの比較をすれば相互作用予測モデルの構築が可能と考えられる．③の原子単位のデータについては，各元素のNN出力を原子特性値とみなして，直接FMO法で得られた教師データと比較をすれば原子特性予測モデルが構築できると考えられる．したがって，各元素のNNの出力を「何とみなす」のか，さらに「みなしたものに応じて総和を取る原子の範囲を調整する」ことで，柔軟に様々な予測モデルを構築することが可能なフレームワークであると言える．

3 原子電荷予測モデルの構築

ここでは，筆者らのグループが報告した原子電荷予測モデル [17]を紹介する．FMOデータを教師データとした原子電荷予測モデルを用いることで，周辺環境に応じた分極・電荷移動を高速に予測可能になると期待される．したがって，時々刻々と変化するタンパク質構造の原子電荷をon the flyで周辺環境に応じて高速に更新するような応用が期待され，静電相互作用の高精度化に寄与すると考えられる．原子電荷予測モデルの教師データについては，FMO創薬コンソーシアムFMODD(FMO Drug Design Consortium) [23]と連携することで，Gaussian等の汎用ソフトでは作成が困難なサイズのタンパク質も含めたデータを作成・利用した．具体的には，polyQ10，TrpCage，BRD2-BD2(BRD2のBD2ドメイン，以後BRD2)について，Amberを用いたMD計算を実施し，それぞれ10,000，10,000，1,000スナップショット構造を作成した．各対象の原子数とFMO計算コストを鑑み，BRD2のみ1,000構造とした．各スナップ書構造に対して，ABINIT-MP [24, 25]を用いたFMO計算(FMO2-HF/6-31G*)によりRESP電荷を算出し教師データとした．MD計算で得られた構造をそのままFMO計算したため一部の構造でFMO計算が収束しなかったものの，FMO計算が収束した各構造のデータをTrain，Validation，Testに分割して学習に用いた．なお，データ分割に際しては，各構造のRMSDを用いたk-meansクラスタリングを行い，得られたクラスタの代表構造及び追加で構造サンプリングを行った構造をTestとした．残りの構造を8:2に分割してTrain，Validationとした(Table 1)．さらに，学習済みモデルの精度を検証するための追加テストデータセットとして，polyQ10とTrp-Cageについては1000 psごとに，BRD2については50 psから500 psごとに追加サンプリングを行った． 最後に，溶媒分子を除去したサンプル分子構造を用いてFMO計算を行った．FMO計算の結果，polyQ10では1000サンプル，TrpCageは989サンプル，BRD2は180サンプルが収束し，追加テストデータとして採用された．

Table 1 Details of the dataset used to build the atomic charge prediction model. Note that the numbers in the table do not include water molecules.

	polyQ10	TrpCage	BRD2
PDBID	2OTU	1L2Y	5IBN
Residues (Atoms)	10(173)	20(304)	111(1,811)
MD snapshots	10,000	10,000	1,000
FMO data	10,000	9,805	902
Train	7,997	7,840	718
Validation	2,000	1,961	180
Test	3+1,000	4+989	4+180

学習については，いずれの系においてもValidationデータに対するRMSEが0.10程度となるまで行った．学習済みの原子電荷予測モデルを用いてクラスタ代表構造のTestデータに対する予測を行った結果をTable 2に示す．また，最もR²の高かったTestデータに対する真値‐予測値プロットはFigure 1のとおりである．原子数の少ないpolyQ10が最も予測精度が高くなっているが，1,811原子から成るBRD2においても，R²で0.9に迫る精度を達成しており，高い予測精度のモデル構築に成功した．ただし，系全体の電荷の総和であるNet chargeを見ると，BRD2についてはNet chargeのブレが大きいため，例えば力場の静電相互作用計算の際に用いる場合には，Net chargeが正しくなるように全ての原子に等しく一定値を加えるような補正を加える等の手立てが必要と考えられる．

Table 2 The coefficient of determination (R²) and the net charge for test data at the cluster center for each target are presented. The net charge represents the sum of the atomic charges of all atoms. The correct values for polyQ10, TrpCage, and BRD2 are 0, 1, and 1, respectively.

	polyQ10		TrpCage		BRD2
Test	R2	Net charge	R2	Net charge	R2	Net charge
1	0.982	-0.31	0.928	1.50	0.896	1.35
2	0.960	-0.23	0.926	0.97	0.888	3.63
3	0.955	0.46	0.917	0.51	0.886	-0.48
4			0.900	1.38	0.873	4.60

Figure 1

 Comparison of FMO calculated and NN predicted values for RESP charges for the structure with the highest coefficient of determination (R²) in the Test data.

4 相互作用予測モデルの構築

次に，相互作用(IFIE)予測モデルの現状 [26]を紹介する．相互作用予測モデルは，ドッキングポーズの評価，動的FMO解析の高速化，IFIE向上を目的変数とした分子生成モデルなど，多岐にわたる応用が期待される．筆者らのグループでは，CDK2を対象として相互作用モデルの基礎検討を進めている．教師データはLINC(Life Intelligence Consortium) [27]とFMODD [23]の連携のもとで生成された．鋳型となるCDK2の構造を選定し，CDK2との複合体構造がpdbに登録されているリガンド(H, C, N, O, Sのみで構成)を80種選定，各リガンドに対してrDockで結合ポーズを20種類生成した後，リガンド周辺を切り出してFMO計算を行うという手順をとった．切り出し方については，FMO計算のコストを鑑みてリガンド周辺9Å以内のアミノ酸残基を対象とする場合に限定して検討した．FMO計算はABINIT-MP [24, 25]を用いてFMO2-MP2/6-31G*レベルで実施した．80種のリガンドのうち，65種をTrain，Validationとし，15種をTestとした．ここでは，リガンドに対する各アミノ酸残基のIFIEの総和(IFIE-Sum)を予測するモデルの構築を行った．

15種のTestリガンドについて，各20ポーズのIFIE-Sumを予測した結果をFigure 2に示す．決定係数は0.60と未知のTestリガンドの様々な結合ポーズに対して良好な予測精度を達成した．特筆すべきはその計算速度であり，本モデルを用いることで一般のノートPCでも1秒とかからずに結果を得ることができる．今後はより実用的なモデルにブラッシュアップすべく，更なる精度向上および対応元素の増強を進める予定である．

Figure 2

 Comparison of FMO calculated and NN predicted values of IFIE-Sum for 20 poses of each of 15 Test ligands of CDK2 (total 300 data).

5 機械学習力場の構築

最後に，機械学習力場構築の現状を紹介する．機械学習力場については筆者らのグループで作成した独自スクリプトのベンチマークを水のみの系で行い，その後，小タンパク質やRNAの検討へと展開している．水のみの系でのベンチマーク結果はFigure 3に示すとおりであり，既報 [28]と同程度の精度を達成している．水のみの系での検証においては，Qunatum-ESPRESSOという周期系の平面波基底を用いた第一原理計算プログラムを用いてデータ生成を行った．

Figure 3

 Benchmark of machine learning forcefield training scripts for water-only system.

現在は，FMO法によりTrpCage，環状ペプチド，RNAといった小~中規模の生体高分子・生体関連分子のデータ生成を進めている．いずれの系においてもHDNNPとACSFsの組み合わせにより学習可能であることは検証済みであり，精度向上に向けた取り組みを推進しているところである．

6 結言

本稿では，筆者らのグループで開発しているFMOデータを用いた機械学習モデルについて概説した．FMO法は原子単位・フラグメント単位・タンパク質単位という3階層にわたる量子化学計算データを生成できる稀有な手法である．それぞれの階層のデータに対する機械学習モデルは，それぞれに有意義な応用が期待されるものである．今後の更なる展開に注目されたい．また，筆者のグループ以外でも，MDトラジェクトリのデータを用いてIFIEを予測するモデルの構築 [29]や，FMO計算で得た分子間相互作用を予測するモデルを用いた高分子相溶性を表すFlory-Hugginsのχパラメータ算定 [30]も報告されている．今後もFMOデータを用いた様々な機械学習モデルの開発が一層盛んになっていくことが期待される．

謝辞

本研究の一部は，AMEDの創薬等先端技術支援基盤プラットフォーム(BINDS)(JP24ama121030)の支援を受けた．FMO計算にはスーパーコンピュータ「富岳」(課題番号: hp240114)および理研所有のスーパーコンピュータ「HOKUSAI」を用いた．また，本研究は，LINCのプロジェクトの一部として実施された．また，FMO計算については，FMO創薬コンソーシアムFMODDとLINCとの連携により実施された．LINC及びFMODDの関係者の皆様にこの場を借りて御礼申し上げる．

参考文献

• [1]   W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, et al., J. Am. Chem. Soc., 117, 5179 (1995). doi:10.1021/ja00124a002
• [2]   W. L. Jorgensen, J. Tirado-Rives, J. Am. Chem. Soc., 110, 1657 (1988). doi:10.1021/ja00214a001 PMID:27557051
• [3]   A. D. MacKerell, Jr., N. Banavali, N. Foloppe, Biopolymers, 56, 257 (2000-2001). PMID:11754339
• [4]   C. Oostenbrink, A. Villa, A. E. Mark, W. F. Van Gunsteren, J. Comput. Chem., 25, 1656 (2004). doi:10.1002/jcc.20090 PMID:15264259
• [5]   J. Behler, M. Parrinello, Phys. Rev. Lett., 98, 146401 (2007). doi:10.1103/PhysRevLett.98.146401 PMID:17501293
• [6]   S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, et al., Nat. Commun., 13, 2453 (2022). doi:10.1038/s41467-022-29939-5 PMID:35508450
• [7]   S. Takamoto, S. Izumi, J. Li, Comput. Mater. Sci., 207, 111280 (2022). doi:10.1016/j.commatsci.2022.111280
• [8]   A. P. Bartók, M. C. Payne, R. Kondor, G. Csányi, Phys. Rev. Lett., 104, 136403 (2010). doi:10.1103/PhysRevLett.104.136403 PMID:20481899
• [9]   G. Kresse, J. Furthmüller, Phys. Rev. B Condens. Matter, 54, 11169 (1996). doi:10.1103/PhysRevB.54.11169 PMID:9984901
• [10]   G. Kresse, J. Furthmüller, Comput. Mater. Sci., 6, 15 (1996). doi:10.1016/0927-0256(96)00008-0
• [11]   G. Kresse, J. Hafner, Phys. Rev. B Condens. Matter, 47, 558 (1993). doi:10.1103/PhysRevB.47.558 PMID:10004490
• [12]   P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, et al., J. Phys. Condens. Matter, 29, 465901 (2017). doi:10.1088/1361-648X/aa8f79 PMID:29064822
• [13]   P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, et al., J. Phys. Condens. Matter, 21, 395502 (2009). doi:10.1088/0953-8984/21/39/395502 PMID:21832390
• [14]Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; et al. Gaussian 16 Rev. C.01. 2016.
• [15]   NequIP; https://github.com/mir-group/nequip.
• [16]   n2p2; https://compphysvienna.github.io/n2p2/index.html.
• [17]K. Kato, T. Masuda, C. Watanabe, N. Miyagawa, H. Mizouchi, S. Nagase, J. Chem Inf. Model. 2020, 60 (7), 3361-3368, Article. doi:10.1021/acs.jcim.0c00273
• [18]   O. T. Unke, M. Stöhr, S. Ganscha, T. Unterthiner, H. Maennel, S. Kashubin, et al., Sci. Adv., 10, eadn4397 (2024). doi:10.1126/sciadv.adn4397 PMID:38579003
• [19]   K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett., 313, 701 (1999). doi:10.1016/S0009-2614(99)00874-X
• [20]   J. Behler, J. Chem. Phys., 134, 074106 (2011). doi:10.1063/1.3553717 PMID:21341827
• [21]   M. Gastegger, L. Schwiedrzik, M. Bittermann, F. Berzsenyi, P. Marquetand, J. Chem. Phys., 148, 241709 (2018). doi:10.1063/1.5019667 PMID:29960372
• [22]   N. Artrith, A. Urban, G. Ceder, Phys. Rev. B, 96, 014112 (2017). doi:10.1103/PhysRevB.96.014112
• [23]   FMO Drug Design Consortium. https://fmodd.jp/
• [24]   Y. Mochizuki, T. Nakano, K. Sakakura, Y. Okiyama, H. Watanabe, K. Kato, et al., The ABINIT-MP Program. In Recent Advances of the Fragment Molecular Orbital Method, SpringerLink, 2021; pp 53-67.
• [25]S. Tanaka, Y. Mochizuki, Y. Komeiji, Y. Okiyama, K. Fukuzawa, Phys. Chem. Chem. Phys., 2014, 16 (22), 10310-10344. doi:10.1039/C4CP00316K
• [26]   R. Kita, C. Watanabe, M. Ohta, N. Tanimura, K. Okuwaki, M. Ikeguchi, et al., Proceedings of the Annual Conference of JSAI 2024, JSAI2024, 4Q3IS2d05-04Q03IS02d05. doi:10.11517/pjsai.JSAI2024.0_4Q3IS2d05
• [27]   Life Intelligence Consortium. https://linc-ai.jp/
• [28]   T. Morawietz, A. Singraber, C. Dellago, J. Behler, Proc. Natl. Acad. Sci. USA, 113, 8368 (2016). doi:10.1073/pnas.1602375113 PMID:27402761
• [29]   S. Tokutomi, K. Shimamura, K. Fukuzawa, S. Tanaka, Chem. Phys. Lett., 757, 137883 (2020). doi:10.1016/j.cplett.2020.137883
• [30]   H. Doi, S. Matsuoka, K. Okuwaki, R. Hatada, S. Minami, R. Suhara, et al., Jpn. J. Appl. Phys., 62, 070901 (2023). doi:10.35848/1347-4065/ace575