Learning Organotransition Metal Reactions Using Graph Neural Networks

Abstract

深層学習による反応予測は，人が識別パターンを設計する必要がないため最近注目されている．しかし有機遷移金属反応は一見複雑で，深層学習の適応例は殆どなかった．我々は，論文や特許などから集めた数万件の実験反応を，素反応に分解したデータベースを構築し，グラフニューラルネットワークを用いて学習させた．最高で97.1%の精度で反応が予測できた．

Translated Abstract

深層学習による反応予測は，人が識別パターンを設計する必要がないため最近注目されている．しかし有機遷移金属反応は一見複雑で，深層学習の適応例は殆どなかった．我々は，論文や特許などから集めた数万件の実験反応を，素反応に分解したデータベースを構築し，グラフニューラルネットワークを用いて学習させた．最高で97.1%の精度で反応が予測できた．

1 はじめに

近年様々な分野に深層学習が応用され，化学でも物性予測，望みの性質を持つ化合物の構造生成，有機化学反応の生成物予測が報告されている [1]．有機遷移金属は産業界でも触媒として広く使われ，全数把握が困難な程，非常に多くの実験，理論研究が報告され，新規反応研究も活発に続いている．しかし深層学習は有機遷移金属反応へ殆ど応用されていない．有機遷移金属は環を多数含む事が多く，構造を文字列で表すSMILES記法が複雑になる．HCNOFなどの典型元素と異なり遷移金属では複数の原子価を取るため，SMILESを扱うライブラリーで問題が起きがちである．そのため，深層学習に必要な良い反応データベースがほとんど無いことが原因である．

深層学習による化学反応の予測では，分子のLewis構造式をグラフとみなしグラフニューラルネットワーク(GNN)を用いるアプローチ [2]と，分子構造式をSMILES記法の文字列で表現し，自然言語処理モデルを用いるアプローチ [3]が発展してきた．後者のアプローチでは，言語のニューラル翻訳モデルを利用するためモデル開発が容易であり，データベースの前処理も簡単である．しかし3次元物体である分子を1次元文字列で表すのは不自然である．また追加の条件(反応条件など)や補助知識(官能基の反応性)をモデルに加えるのは難しい．一方前者のアプローチでは，グラフで分子構造を直接表現するため，ネットワークと実際の分子の対応が明らかで，有機遷移金属反応によりふさわしい．そこで本研究では，有機遷移金属反応のデータベースとGNNを用いた機械学習モデルを構築し，その能力を調べることにした．

2 方法

データベース:学習に用いる実験反応データは，米国特許情報から自動収集された化学反応 [4]のうち，1 個以上の遷移金属原子を含むもの約5.5万件 (Patent data)と，CAS (CAS, a division of the American Chemical Society)から提供して頂いた，重要論文が報告した反応約8千件 (Journal data)である．これをReaction Decoder Tool [5]を用いて反応物と生成物の原子のマッピングを行った．有機遷移金属の触媒反応では，何種類かの基本反応(素反応)を順に行い生成物ができる．素反応はより単純であり，ある種の素反応は複数回現れ，学習に使える反応数も増えるため，素反応を学習させることにした．有機遷移金属化学の教科書 [6]を参考に実験反応を素反応に分解し，学習に適したデータベースを構築した．

モデル:分子をLewis構造式のグラフで表した．そのnodeには元素を，edgeには結合種類を表すone-hotラベルを付けた．SMILES記法の反応データをこれら分子グラフに変換した．本研究のグラフニューラルネットワーク(GNN)は，入力された反応物に対し，結合ができるか切れる原子対を予測することで，生成物を予測する(Figure 1)．まず，GNNの畳み込み層では近傍の原子N(i)の特徴ベクトルを，結合情報を加味しながら集約して，原子iの特徴ベクトルを更新していく[Figure 1, 式(1), (2)]．

Figure 1.

 Proposed neural network model.

ここで，Θ_kは学習可能なパラメータを示している．   

${a_i}^{\text{'}}={\Theta }_1{a}_i+{\sum }_{j\in N\left(i\right)}h\left(b_{i,j}\right)a_j$(1)

  

$h\left(b_{i,j}\right)={\Theta }_3\left(ReLU\left({\Theta }_2{b}_{i,j}\right)\right)$(2)

この特徴ベクトルから，各結合が切断される確率や，各原子対が新たに結合する確率を計算する．反応性の予測には，HSAB則とフロンティア軌道理論に着目した．Drago と WaylandはHSAB則を使い，酸塩基相互作用のエンタルピーを式(3)で近似できることを報告した [7]．この式はベクトル$a_A=\left(E_A, C_A\right)$と$a_B=(E_B, C_B)$の内積である．フロンティア軌道理論によれば，ドナー分子Dとアクセプター分子Aの相互作用強度は積分$\int {\phi }_{HOMO}{\phi }_{LUMO}d\tau$で与えられる．ここで${\phi }_{HOMO}$はDのHOMO，${\phi }_{LUMO}$はAのLUMOである．MO ${\phi }_{HOMO}$を原子軌道(AO)${\chi }_k$と対応するMO係数$C_k^{HOMO}$で展開すると積分は$\sum _{k,l}{C}_k^{HOMO}{C}_l^{LUMO}{V}_{k,l}$となる．適切な混成AOを用いれば，この相互作用行列$V_{k,l}$は対角優位になる．したがって，HOMO-LUMO相互作用は原子対の寄与の和として近似される[式(4)]．   

$-∆H_{AB}=E_A{E}_B+C_A{C}_B$(3)

  

$\sum _{k,l}\left(C_k^{HOMO},C_k^{LUMO}\right)\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}{C}_l^{HOMO}\\ C_l^{LUMO}\end{array}\right)=\sum _{k,l}{a}_k^{}\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)a_l^T$(4)

  

$a_k=\frac{1}{\sqrt{2}}\left(C_k^{HOMO}+C_k^{LUMO},C_k^{HOMO}-C_k^{LUMO}\right)$(5)

そこで，原子の特徴ベクトル$a_i$の内積が，各結合の反応性を表すと仮定した．原子A，Bが金属Mから脱離してA，B間に結合ができる反応性$f^{eli}(x_A, x_B)$を式(6)，原子A，B間の結合を切断して金属Mが付加する反応性$f^{add}(x_A, x_B)$を式(7)のように仮定した．ここでΘ_Mはバイアス項を持つ線形層を示す．この層には遷移金属原子の特徴を入力し，高配位金属原子の付加反応性が下がることを表現させた．   

$f^{eli}(x_A, x_B) = x_A\bullet x_B-x_M\bullet x_A- x_M\bullet x_B$(6)

  

$f^{add}(x_A, x_B) = - f^{eli}(x_A, x_B) - {\Theta }_M{x}_M$(7)

この式が予測した最も反応性の高い原子対が実験値を再現するように，交差エントロピーを最小化してGNNの重みを学習した．機械学習にはPyTorchを用いた．

3 結果と考察

Journal data及びPatent dataのtest setでの予測精度はそれぞれ95.4%及び97.1%となった(Table 1)．また，反応別では特徴的な反応剤を用いるものは高精度だったが，C-H結合を変換する反応の精度はより低かった．C-H結合はありふれているが，反応するC-H結合は珍しく，反応性が低く見積もられている可能性がある．基質の金属への付加と金属からの脱離の2パターンを比較すると，脱離のほうが高精度だった．これは金属を手掛かりに反応位置を絞り込めるためと説明できる．同じ基質に対して異なる金属触媒を用いた時，金属種によって予想部位が異なった(Figure 2)．つまりモデルは金属の特徴を学習したことが示唆される．また，ノードの特徴ベクトルを主成分分析したところ，原子の種類や分子の部分構造に応じてクラスターを形成していた．これらを認識して反応性を計算していると考えられる．予測に失敗した例では，反応活性部位が複数含まれる分子で，正解以外の活性部位を予測した例が多かった．本モデルでは反応性を制御する配位子を無視しているため，予測が難しい例があることを端的に表している．

Table 1. Prediction accuracy of selected reactions.

Figure 2.

 Visualization of bond reactivity.

4 結論

本研究では，素反応に分解した有機遷移金属反応のデータベースおよびGNNを用いたモデルを構築し，化学反応予測を行った．全体の精度は最大97.1%となり，特徴的な反応剤を用いる反応を非常に高精度で予測できることが分かった．モデルの振る舞いや原子の特徴ベクトルの主成分分析から，原子の種類や分子の部分構造を理解して反応性を計算していることが示唆された．

謝辞

CASのデータ収集にご尽力下さった中村様，Dr. Farr, Mr. Thompson, Dr. Sonに厚く御礼申し上げます．この研究は科研費(7321K12027)と名古屋大学融合フロンティアフェローシップ(JPMJFS2120)の補助を受けています．

参考文献

• [1]   M. Meuwly, Chem. Rev., 121, 10218 (2021). doi:10.1021/acs.chemrev.1c00033 PMID:34097378
• [2]   V. R. Somnath, C. Bunne, C. W. Coley, A. Krause, R. Barzilay, ArXiv, vol. abs/2006.07038, 2020.
• [3]   P. Schwaller, T. Laino, T. Gaudin, P. Bolgar, C. A. Hunter, C. Bekas, et al., ACS Cent. Sci., 5, 1572 (2019). doi:10.1021/acscentsci.9b00576 PMID:31572784
• [4]   D. Lowe, Chemical Reactions from US Patents (1976-Sep2016), Figureshare dataset, (accessed 2023–1-16).doi:10.6084/m9.figshare.5104873.v1
• [5]   S. A. Rahman, G. Torrance, L. Baldacci, S. Martínez Cuesta, F. Fenninger, N. Gopal, et al., Bioinformatics, 32, 2065 (2016). doi:10.1093/bioinformatics/btw096 PMID:27153692
• [6]   O. M. Chemistry, From Bonding to Catalysis, J. F. Hartwig (University Science Books, 2010).
• [7]   R. S. Drago, B. B. Wayland, J. Am. Chem. Soc., 87, 3571 (1965). doi:10.1021/ja01094a008