2021 Volume 20 Issue 3 Pages 106-111
We have developed a novel framework for inferring a chemical graph. Given an ANN that
maps a feature vector
構造活性相関,構造物性相関に人工ニューラルネットワーク(ANN)などの機械学習を用いた場合に,指定した物性値を予測値として呈する分子構造を求める逆解析問題を数理計画法の利用により数学的に解く新規手法を開発した.本稿ではこの方法の最新版 [1]の概要について紹介し,この文献で扱われていない新たな物性を用いた計算機実験結果について報告を行う.
提案法は,物性値を予測する関数を機械学習で構築する段階1-3,ならびに,物性値予測の逆問題を解く段階 4, 5 から成り,後半の逆解析は数理計画法の手法を用い数学的に正確に解決する.また,目標の化学グラフの形状を指定・制限するために構造仕様と呼ぶ記述ルールを導入する.
2.1 化学グラフ本稿でグラフは連結単純無向グラフを指す.次数(隣接点の個数)が1である点を葉と呼ぶ.閉路を持たないグラフを木と呼び,枝が直列につながったグラフをパスと呼ぶ.化合物は化学グラフ(点にラベル,枝に結合度が付与された無向グラフ)として扱う(立体構造の情報は含まない).非環式(木状構造)化学グラフは1点が根として指定されているとき化学根付き木と呼ぶ.
2.2 機械学習段階1: 対象とする性質πの実測値
段階2: 特徴関数
段階3: データセット
化学グラフの構造を外部と内部の二つに分けて捉える二層モデルを提案する.
外部と内部の境目を定める整数値ρを分岐高さ(branch-height)と呼び, 通常ρ = 2と設定している.水素の表示を無視した化学グラフ
(a) A chemical graph
内部の特徴化 内部の2点
外部の特徴化 データセット
特徴関数は枝構成と外縁構成を主要な記述子とするベクトルと定める.
2.4 構造仕様逆解析で求めたい対象の化学グラフの形状を限定する仕方について考える.以下のように,(i)-(iii)で内部を,(iv)で外部を定める.
(i) 二層モデルでの内部の大まかな形を表現するグラフ
(ii) 種グラフの枝ごとに,削除あるいはパス(これを純粋パスと呼ぶ)への置き換えを行う.
(iii) 次に,いくつかの点からパス(これを葉パスと呼ぶ)を外側へ伸ばすことで,グラフを拡大させる.こうして得られたグラフを目標の化学グラフの内部とする.
(iv) 最後に,内部の点
Figure 1 (b)に種グラフ
構造仕様の設定の仕方により化合物空間における探索範囲の広さを自在に制御することができる.
2.5 逆解析の仕組み特徴関数
段階4:
特徴関数f,予測関数η,構造仕様σを模倣するMILPを定式化しておく. 所望の物性値
段階5: 段階4で構築した化学グラフG†を元にして,内部の枝の配列,外部の外縁木の配列を並び変える操作により,同じ目標値y*,特徴関数,構造仕様を 満足する構造異性体G*を動的計画法アルゴリズムにより多数(あるいは指定の個数まで)生成する.
以下の計算機実験ではパソコン(3.0GHz Core i7-9700, Memory:16GB RAM DDR4)を使用し,ANNの構築にはscikit-learn (0.23.2)を用い,MILPのソルバーとしてはCPLEX (12.10)を利用した.
段階3:文献 [1]では実験されていない以下の 七種の物性に対し,化学グラフ(立体情報は含まない)及び物性の実測値 を用いて段階3の実験を行った. 蒸気密度 (vapor density) [2],旋光 (optical rotation) [2],絶対零度での内部エネルギー (internal energy at 0K) [2],比熱容量 (heat capacity at 298.15K) [3],等方性分極率 (isotropic polarizability) [3],等圧熱容量(液体・固体) (isobaric heat capacities liquid/solid) [4].各物性ごとに,予備実験として,記述子の個数を絞り込み,ANNのアーキテクチャを複数個の候補から一つ選定しておいた.性能試験では10セットの交差検定を行った(1セットの検定ではデータセットD を4対1の比で教師データ,テストデータとして5回使う).この結果,七種の物性に対してそれぞれの テストデータの決定係数R2の中央値は0.90∼0.99であった.
段階4:段階3で得たANNを用いて指定の構造仕様を満足する化学グラフG†をMILPの解として一つ求める.実験用に次のような例題I1, I2, I3 を用意した(仕様の詳細は [1]を参照; CIDはPubChemの化合物id), 例題I1: 二つの環が1点を共有した構造に化学根付き木が付随した化合物(種グラフはFigure 2 (a)を参照); 例題I2: 化合物CID24822711の環構造部分(Figure 2 (c)の灰色部分)を内部に有し,外縁木には化合物CID59170444 (Figure 2 (d))の外縁木のみを用いた化合物; 例題I3: 二つの単環式化合物CID10076784,CID44340250 (Figure 2 (e), (f)参照) の特徴関数の中間の特徴ベクトルを呈する単環式化合物.
(a) A seed graph
指定した等圧熱容量(液体) (IhcLiq) の目標値y*に対して 各例題のMILPを解いた結果として,MILPの変数の個数,制約式の本数,MILPの求解時間(秒),得られた化学グラフG†の非水素原子数を表1に示す.非水素原子数が50前後でも10秒以下で例題を解くことができている.例題I1, I2, I3に対する逆解析の解である化学グラフG†を それぞれFigure 2 (b), (g), (h)に図示する.例題I1に対するG†は単環の化学グラフである.例題I2に対するG†は,内部に化合物CID24822711の環構造部分を有し, 外部にはCID59170444の外縁木のみを有している.I3に対するG†は CID10076784とCID44340250の間にある構造を有している.いずれの解の形状も構造仕様の指定に従っている.このように構造仕様のルールを工夫して利用することで逆解析において構築を望む化合物の形状を計画的に設定することができる.
複数の物性の目標値を同時に満たす化学グラフG† を推定する設定も容易に扱える.
段階5:同じ目標値y*,構造仕様を満足するG†の構造異性体G*を生成する. 提案アルゴリズムは答えとなるG*を列挙せずに総数の下界も計算できる.実験では高々100個まで生成させている.アルゴリズムの計算時間(秒),同じ仕様を満たすG†の構造異性体G*の総数の下界,100個を上限として実際に生成した異性体の個数をTable 2に表す.例題 I3の計算では106を超える個数の構造異性体の存在を検出できている.G†の構造によっては,枝や外縁木の並び変えに自由度がない場合,条件を満たす別の異性体G*がほとんど見つからないこともある. このようなときは構造仕様あるいは段階4のMILPの条件に若干手を加えて別の解を構築させる方法も設計されている [5].
例題 | 目標値y* | 変数,制約式 の個数 | 計算時間 (秒) | 原子数 |
I1 | 750 | 11054, 8656 | 8.60 | 50 |
I2 | 600 | 7776, 8393 | 7.25 | 42 |
I3 | 1187 | 5304, 6662 | 3.55 | 42 |
例題 | 計算時間 (秒) | 異性体数の下界 | 生成した個数 |
I1 | 0.44 | 60 | 60 |
I2 | 0.12 | 35 | 35 |
I3 | 1.15 | 2.5×107 | 100 |
ANNにより構造と物性の相関が学習できたとき, 任意の形状の分子構造の化合物を推定する方法を提案した.このために以下の独自の概念,技術を開発した. ANNの計算を模倣するMILPの定式化,二層モデルに基づく特徴関数,推定したい分子構造の形状をルールで規定する構造仕様,同じ特徴ベクトル・構造仕様を満足する異性体を列挙する動的計画法アルゴリズム.今後は,線形回帰,決定木,グラフ畳み込みなど他の機械学習に基づく逆解析法を設計していく予定である.