Journal of Computer Chemistry, Japan
Online ISSN : 1347-3824
Print ISSN : 1347-1767
ISSN-L : 1347-1767
Letters (Selected Papers)
[title in Japanese]
Hiroshi NAGAMOCHIJianshen ZHUNaveed Ahmed AZAMKazuya HARAGUCHILiang ZHAOTatsuya AKUTSU
Author information
JOURNAL FREE ACCESS FULL-TEXT HTML

2021 Volume 20 Issue 3 Pages 106-111

Details
Abstract

We have developed a novel framework for inferring a chemical graph. Given an ANN that maps a feature vector x = f ( G ) of a chemical graph G to a predicted value η ( x ) of a chemical property to G, we formulate an integer program that simulates the functions f and η. Given a target value y*, we can infer a chemical graph G such that η ( f ( G ) ) = y * by solving the integer program.

1 はじめに

構造活性相関,構造物性相関に人工ニューラルネットワーク(ANN)などの機械学習を用いた場合に,指定した物性値を予測値として呈する分子構造を求める逆解析問題を数理計画法の利用により数学的に解く新規手法を開発した.本稿ではこの方法の最新版 [1]の概要について紹介し,この文献で扱われていない新たな物性を用いた計算機実験結果について報告を行う.

2 方法

提案法は,物性値を予測する関数を機械学習で構築する段階1-3,ならびに,物性値予測の逆問題を解く段階 4, 5 から成り,後半の逆解析は数理計画法の手法を用い数学的に正確に解決する.また,目標の化学グラフの形状を指定・制限するために構造仕様と呼ぶ記述ルールを導入する.

2.1 化学グラフ

本稿でグラフは連結単純無向グラフを指す.次数(隣接点の個数)が1である点を葉と呼ぶ.閉路を持たないグラフを木と呼び,枝が直列につながったグラフをパスと呼ぶ.化合物は化学グラフ(点にラベル,枝に結合度が付与された無向グラフ)として扱う(立体構造の情報は含まない).非環式(木状構造)化学グラフは1点が根として指定されているとき化学根付き木と呼ぶ.

2.2 機械学習

段階1: 対象とする性質πの実測値 a ( G ) が既知である化学グラフ G を集めたデータセットDを用意する.

段階2: 特徴関数 f として化学グラフ G の特徴を抽出したいくつかの 記述子からなる特徴ベクトル f ( G ) を定義する. Kを使用する記述子の個数としたとき, f は化合物空間から特徴空間 K への関数である.

段階3: データセット D を教師データとして,ANN,決定木などの機械学習を用いて, 特徴ベクトル f ( G ) と実測値 a ( G ) の関係を学習させ,特徴空間 K から実数値空間 への予測関数 η を構築する.これにより新規化学グラフGに対する物性値 a ( G ) の予測値として η ( f ( G ) ) が利用できる.

2.3 二層モデルに基づく特徴関数

化学グラフの構造を外部と内部の二つに分けて捉える二層モデルを提案する. 外部と内部の境目を定める整数値ρを分岐高さ(branch-height)と呼び, 通常ρ = 2と設定している.水素の表示を無視した化学グラフ において, 葉から高さρまでの化学根付き木を外部と呼び,残りを内部と呼ぶ.Figure 1 (a)の化学グラフの例に対し 外部は灰色部の8個の化学根付き木であり,それぞれは点 を根とする木 G とみなし,これらをGの外縁木(fringe-tree)と呼ぶ. T u i から外縁木の根以外の点を除いた残りの部分を G の内部と定める.

Figure 1.

 (a) A chemical graph y * ; (b) A seed graph G ; (c) A set G * of chemical rooted trees; (d) A graph G * ; (e) A graph G .

内部の特徴化 内部の2点 G をつなぐ枝 G の特徴量として,点 u , v の元素a,点 u の隣接点数d,点 u の元素b,点 v の隣接点数h,ならびに 枝 v の結合度 h からなる組 e を用いる(これを枝構成(edge-configuration)と呼ぶ).

外部の特徴化 データセット m 内の化学グラフが持つ外縁木を全種類,事前に調べておきコード化しておく(物性にもよるがρ = 2であれば数十∼200種程度). 外部の特徴量としては, ( a d , b h , m ) の外縁木の頻度をコードで表したものを用いる(これを外縁構成(fringe-configuration)と呼ぶ).

特徴関数は枝構成と外縁構成を主要な記述子とするベクトルと定める.

2.4 構造仕様

逆解析で求めたい対象の化学グラフの形状を限定する仕方について考える.以下のように,(i)-(iii)で内部を,(iv)で外部を定める.

(i) 二層モデルでの内部の大まかな形を表現するグラフ D (種グラフ(seed graph)と呼ぶ)を与える.種グラフの枝は,そのまま利用するもの,削除を許すもの,パスに 取り換えることを許すものが指定されている.

(ii) 種グラフの枝ごとに,削除あるいはパス(これを純粋パスと呼ぶ)への置き換えを行う.

(iii) 次に,いくつかの点からパス(これを葉パスと呼ぶ)を外側へ伸ばすことで,グラフを拡大させる.こうして得られたグラフを目標の化学グラフの内部とする.

(iv) 最後に,内部の点 ρ = 2 に化学根付き木をその点を根とする外縁木Tuとして添加し,目標の化学グラフとする.

Figure 1 (b)に種グラフ G C の例を示す.破線の枝eE(≥2),点線の枝eE(≥1)はそれぞれ長さ2以上,1以上の純粋パスに交換することを許し,灰色の枝eE(0/1)は削除することを許し,黒実線の枝eE(=1)はそのまま使用することを意味している. 種グラフ u の指定の他,(iv)の外縁木として使うことを許す化学根付き木の集合 T u の指定,(ii)の純粋パスや(iii)の葉パスの長さの上下限などの指定をルールとして持たせたものが構造仕様である.Figure 1 (a)の化学グラフ G C Figure 1 (b)の種グラフ Figure 1 (c)の集合 e E ( 1 ) およびその他の条件 (詳しくは [1]参照)から成る構造仕様σを満足する例となっている.まず, の枝 G C が純粋パス F に置き換えられ, Figure 1 (d)に示すグラフ G が得られる.次に,点 G C を根とする葉パス F が付与され, Figure 1 (e)に示すグラフ が内部構造として得られる.最後に, G C の枝に結合度,点に元素を与え, a i , i = 1 , 2 , 3 , 4 , 5 から選ばれた化学木を外縁木 P a i として添加することでFigure 1 (a)の化学グラフ G 1 が得られる.

構造仕様の設定の仕方により化合物空間における探索範囲の広さを自在に制御することができる.

2.5 逆解析の仕組み

特徴関数 u i , i = 5 , 18 , 22 は元素や隣接元素対の含有頻度など簡単な記述子からなるものと仮定する.機械学習にANNや決定木を用いる場合,特徴関数 Q u i や予測関数 G 2 の逆関数の計算を正確に行うアプローチとして 混合整数線形計画問題 (Mixed Integer Linear Programming,MILP) を用いることができる.化学グラフGから F を計算する過程,ならびに,ベクトルx ∈ Kから T u i , i = 12 , 14 , 15 , 21 , 23 , 24 , 27 , 28 を計算する過程をMILPとして模倣(実数・整数変数の線形不等式の集まりとして表現)することができる. 逆解析においては, G を司る変数 f の値を指定の物性値 f に 固定したMILPを解くことで, η を満たす特徴空間のベクトル G および, f(G) =x*なる化学グラフ f ( G ) が一つ得られる.構造仕様σもまた MILPとして模倣できるように導入してある.提案する逆解析法は以下のようになる.

段階4: 特徴関数f,予測関数η,構造仕様σを模倣するMILPを定式化しておく. 所望の物性値 η ( x ) に対してこのMILPを解き,および 構造仕様を満たす化学グラフGを一つ求める.

段階5: 段階4で構築した化学グラフGを元にして,内部の枝の配列,外部の外縁木の配列を並び変える操作により,同じ目標値y*,特徴関数,構造仕様を 満足する構造異性体G*を動的計画法アルゴリズムにより多数(あるいは指定の個数まで)生成する.

3 結果

以下の計算機実験ではパソコン(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)参照) の特徴関数の中間の特徴ベクトルを呈する単環式化合物.

Figure 2.

 (a) A seed graph G * ; (b) A chemical graph 10 6 for I1; (c) CID 24822711; (d) CID 59170444; (e) CID 10076784; (f) CID 44340250; (g) A chemical graph G * for I2; (h) A chemical graph / / for I3.

指定した等圧熱容量(液体) (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].

Table 1.  Results of Stage 4 for property IhcLiq.
例題 目標値y* 変数,制約式 の個数 計算時間 (秒) 原子数
I1 750 11054, 8656 8.60 50
I2 600 7776, 8393 7.25 42
I3 1187 5304, 6662 3.55 42
Table 2. Results of Stage 5 for property IhcLiq.
例題 計算時間 (秒) 異性体数の下界 生成した個数
I1 0.44 60 60
I2 0.12 35 35
I3 1.15 2.5×107 100

4 まとめ

ANNにより構造と物性の相関が学習できたとき, 任意の形状の分子構造の化合物を推定する方法を提案した.このために以下の独自の概念,技術を開発した. ANNの計算を模倣するMILPの定式化,二層モデルに基づく特徴関数,推定したい分子構造の形状をルールで規定する構造仕様,同じ特徴ベクトル・構造仕様を満足する異性体を列挙する動的計画法アルゴリズム.今後は,線形回帰,決定木,グラフ畳み込みなど他の機械学習に基づく逆解析法を設計していく予定である.

参考文献
 
© 2021 Society of Computer Chemistry, Japan
feedback
Top