An Inverse QSAR Method Based on a Two-Layered Model and Integer Programming

Abstract

A novel framework for inverse quantitative structure–activity relationships (inverse QSAR) has recently been proposed and developed using both artificial neural networks and mixed integer linear programming. However, classes of chemical graphs treated by the framework are limited. In order to deal with an arbitrary graph in the framework, we introduce a new model, called a two-layered model, and develop a corresponding method. In this model, each chemical graph is regarded as two parts: the exterior and the interior. The exterior consists of maximal acyclic induced subgraphs with bounded height, the interior is the connected subgraph obtained by ignoring the exterior, and the feature vector consists of the frequency of adjacent atom pairs in the interior and the frequency of chemical acyclic graphs in the exterior. Our method is more flexible than the existing method in the sense that any type of graphs can be inferred. We compared the proposed method with an existing method using several data sets obtained from PubChem database. The new method could infer more general chemical graphs with up to 50 non-hydrogen atoms. The proposed inverse QSAR method can be applied to the inference of more general chemical graphs than before.

1. Introduction

Computer-aided design of chemical structures is one of the key topics in chemoinformatics. In particular, extensive studies have been done on inverse quantitative structure–activity relationships (inverse QSAR), which seek chemical structures having desired chemical activities under some constraints. In this framework, chemical compounds are usually represented as vectors of real or integer numbers, which are often called descriptors in chemoinformatics and correspond to feature vectors in machine learning. Using these chemical descriptors, various heuristic and statistical methods have been developed for inverse QSAR [1,2,3]. In many of such methods, inference or enumeration of graph structures from a given set of descriptors is a crucial subtask. Although various methods have been developed for that purpose [4,5,6,7], enumeration still remains a challenging task because the number of possible chemical graphs is huge, for example, chemical graphs with up to 30 atoms (vertices) C, N, O, and S, may exceed ${10}^{60}$ [8]. Furthermore, even inference is a challenging task because it is NP-hard (computationally difficult) except for some simple cases [9]. Due to this inherent difficulty, most existing methods for inverse QSAR do not guarantee optimal or exact solutions.

On the other hand, the design of novel graph structures has recently become a hot topic in artificial neural network (ANN) studies, and thus extensive studies have been done for inverse QSAR using ANNs, especially with graph convolutional networks [10]. For example, variational autoencoders [11], recurrent neural networks [12,13], grammar variational autoencoders [14], generative adversarial networks [15], and invertible flow models [16,17] have been applied. Note that QSAR using three-dimensional structures of chemical compounds (3D-QSAR) has also been studied [18]. Particularly, comparative molecular field analysis (CoMFA) has been extensively studied and applied to various molecular design problems [19,20]. In CoMFA, electrostatic potential interaction energies across superimposed molecular structures are used as descriptors and then regression is performed by using the partial least squares (PLS) fitting. Recently, deep neural networks have been applied to 3D-QSAR by combining potential interaction energies with convolutional neural networks [21]. However, in order to apply 3D-QSAR, we need to calculate accurate three-dimensional structures of chemical compounds, which is not a straightforward task.

A novel framework for inferring chemical graphs has recently been developed [22,23] based on ANNs and mixed integer linear programming (MILP), as illustrated in Figure 1. It constructs a prediction function in the first phase and infers a chemical graph in the second phase. The first phase of the framework consists of three stages. In Stage 1, we choose a chemical property $\pi$ and a class $\mathcal{G}$ of graphs, where a property function a is defined so that $a\left(G\right)$ is the value of $\pi$ in $G\in \mathcal{G}$, and collect a data set $D_{\pi }$ of chemical graphs in $\mathcal{G}$ such that $a\left(G\right)$ is available. In Stage 2, we introduce a feature function $f:\mathcal{G}\to {\mathbb{R}}^K$ for a positive integer K. In Stage 3, we construct a prediction function ${\eta }_{\mathcal{N}}$ with an ANN $\mathcal{N}$ that, given a vector $x\in {\mathbb{R}}^K$, returns a value $y={\eta }_{\mathcal{N}}\left(x\right)\in \mathbb{R}$ so that ${\eta }_{\mathcal{N}}\left(f\left(G\right)\right)$ serves as a predicted value to $a\left(G\right)$ for each $G\in D_{\pi }$. Given a target chemical value $y^{*}$, the second phase infers chemical graphs $G^{*}$ with ${\eta }_{\mathcal{N}}\left(f\left(G^{*}\right)\right)=y^{*}$ in the next two stages. In Stage 4, we formulate an MILP that simulates the construction of $f\left(G\right)$ from G and the computation process in the ANN so that given a target value, $y^{*}$, and solve the MILP to infer a chemical graph $G^{†}$ and a feature vector $x^{*}$ such that $f\left(G^{†}\right)=x^{*}$ and ${\eta }_{\mathcal{N}}\left(x^{*}\right)=y^{*}$. In Stage 5, we generate other chemical graphs $G^{*}$ such that ${\eta }_{\mathcal{N}}\left(f\left(G^{*}\right)\right)=y^{*}$ based on the output chemical graph $G^{†}$.

MILP formulations required in Stage 4 have been designed for chemical compounds with cycle index 0 (i.e., acyclic) [23,24], cycle index 1 [25], and cycle index 2 [26]. In particular, Azam et al. [24] introduced a restricted class of acyclic graphs that is characterize by an integer $\rho$, called a “branch-parameter” such that the restricted class still covers most of the acyclic chemical compounds in the database.

Recently, Akutsu and Nagamochi [27] extended the idea to define a restricted class of cyclic graphs, called “$\rho$-lean cyclic graphs”, that covers most of the cyclic chemical compounds in the database. Based on this, they also defined a set of rules for specifying several topological substructures of a target chemical graph in a flexible way in Stage 4 before we solve an MILP. The method has been implemented by Zhu et al. [28], and computational results showed that chemical graphs with around up to 50 non-hydrogen atoms can be inferred. Although the method can infer the class of $\rho$-lean cyclic graphs and specify topological structures of the cyclic part, we still need to introduce a new model to deal with an arbitrary graph and to include a prescribed structure in the acyclic part of a target chemical graph.

In this paper, we introduce a new model, called a two-layered model, for representing the feature of a chemical graph in order to deal with an arbitrary graph in the framework. In the two-layered model, a chemical graph G with a parameter $\rho \ge 1$ is regarded as two parts: the exterior and the interior. The exterior consists of maximal acyclic induced subgraphs with height at most $\rho$ and the interior is the connected subgraph obtained by ignoring the exterior. We define a feature vector $f\left(G\right)$ of a chemical graph G to be the frequency of adjacent atom pairs in the interior and the frequency of chemical acyclic graphs in the exterior. Figure 2 illustrates an example of a chemical graph G. For a branch-parameter $\rho =2$, the interior of the chemical graph G in Figure 2 is obtained by removing the set of vertices with degree 1 $\rho =2$ times, i.e., first remove the set $V_1=\{w_1,w_2,\dots ,w_{14}\}$ of vertices of degree 1 in G, and then remove the set $V_2=\{w_{15},w_{16},\dots ,w_{19}\}$ of vertices of degree 1 in $G-V_1$, where the removed vertices become the exterior-vertices of G and there are eight rooted trees $T_1,T_2,\dots ,T_8$ in the exterior of G.

We also introduce a new set of rules for specifying topological substructures of a target chemical graph G to be inferred so that a prescribed structure can be included in both of the acyclic and cyclic parts of G. The set of rules contains (i) a seed graph $G_{\mathrm{C}}$ as an abstract form of a target chemical graph G; (ii) a set $\mathcal{F}$ of chemical rooted trees as candidates for trees in the exterior of G; and (iii) lower and upper bounds on the number of components in a target chemical graph such as chemical elements, double/triple bounds and the interior-vertices in G. Figure 3a,b illustrates examples of a seed graph $G_{\mathrm{C}}$ and a set $\mathcal{F}$ of chemical rooted trees, respectively. Given a seed graph $G_{\mathrm{C}}$, the interior of a target chemical graph G is constructed from $G_{\mathrm{C}}$ by replacing some edges $a=uv$ with paths $P_a$ between the end-vertices u and v, and by attaching new paths $Q_v$ to some vertices v. For example, the chemical graph G in Figure 2 is constructed from the seed graph $G_{\mathrm{C}}$ in Figure 3a as follows. First replace five edges $a_1=u_1{u}_2,a_2=u_1{u}_3,a_3=u_4{u}_7,a_4=u_{10}{u}_{11}$ and $a_5=u_{11}{u}_{12}$ in $G_{\mathrm{C}}$ with new paths $P_{a_1}=(u_1,u_{13},u_2)$, $P_{a_2}=(u_1,u_{14},u_3)$, $P_{a_3}=(u_4,u_{15},u_{16},u_7)$, $P_{a_4}=(u_{10},u_{17},u_{18},u_{19},u_{11})$ and $P_{a_5}=(u_{11},u_{20},u_{21},u_{22},u_{12})$, respectively, to obtain the subgraph $G_1$ of G that consists of vertices depicted with squares. Next, attach to this graph $G_1$ three new paths, $Q_{u_5}=(u_5,u_{24})$, $Q_{u_{18}}=(u_{18},u_{25},u_{26},u_{27})$, and $Q_{u_{22}}=(u_{22},u_{28})$, to obtain the interior of G in Figure 2. Finally, the chemical graph G in Figure 2 is obtained by attaching eight trees $T_1,T_2,\dots ,T_8$ selected from the set $\mathcal{F}$ and assigning chemical elements and bond-multiplicities in the interior. The frequency of chemical elements and the graph size are controlled with lower and upper bounds on the components in a target chemical graph G. See Section 2.2 for more details on the specification.

We implemented the two-layered model and the results of computational experiments suggest that the proposed method can infer chemical graphs with around up to 50 non-hydrogen atoms.

The paper is organized as follows. Section 2.1 introduces some notions on graphs, a modeling of chemical compounds, and a choice of descriptors. Section 2.2 introduces a method of specifying topological substructures of target chemical graphs in Stage 4. Section 3 reports the results on some computational experiments conducted for chemical properties such as octanol/water partition coefficient, boiling point, melting point, flash point, lipophylicity, and solubility. Section 4 makes some concluding remarks. An MILP formulation used in Stage 4 and a review of the dynamic programming algorithm for generating isomers in Stage 5 can be found in Supplementary Materials. The proposed method/system is available at GitHub https://github.com/ku-dml/mol-infer.

2. Materials and Methods

This section presents mathematical details of our developed method. Readers not interested in mathematical details can skip this section.

2.1. Preliminary

This section introduces some notions and terminology on graphs, a modeling of chemical compounds, and our choice of descriptors.

Let $\mathbb{R}$, $\mathbb{Z}$ and ${\mathbb{Z}}_{+}$ denote the sets of reals, integers and non-negative integers, respectively. For two integers a and b, let $[a,b]$ denote the set of integers i with $a\le i\le b$.

Graphs. Given a graph G, let $V\left(G\right)$ and $E\left(G\right)$ denote the sets of vertices and edges, respectively. For a subset $V^{\prime }\subseteq V\left(G\right)$ (resp., $E^{\prime }\subseteq E\left(G\right))$ of a graph G, let $G-V^{\prime }$ (resp., $G-E^{\prime }$) denote the graph obtained from G by removing the vertices in $V^{\prime }$ (resp., the edges in $E^{\prime }$), where we remove all edges incident to a vertex in $V^{\prime }$ in $G-V^{\prime }$. The rank $\mathrm{r}\left(G\right)$ of a graph G is defined to be the minimum $\left|F\right|$ of an edge subset $F\subseteq E\left(G\right)$ such that $G-F$ contains no cycle. A path with two end-vertices u and v is called a $u,v$-path. An edge $e=u_1{u}_2$ in a connected graph G is called a bridge if the graph $G-e$ obtained from G by removing edge e is not connected, i.e., $G-e$ consists of two connected graphs $G_i$ containing vertex $u_i$, $i=1,2$. For a cyclic graph G, an edge e is called a core-edge if it is in a cycle of G or is a bridge $e=u_1{u}_2$ such that each of the connected graphs $G_i$, $i=1,2$ of $G-e$ contains a cycle. A vertex incident to a core-edge is called a core-vertex of G.

A vertex designated in a graph G is called a root. In this paper, we designated at most two vertices as roots, and denote by $\mathrm{Rt}\left(G\right)$ the set of roots of G. We call a graph G rooted (resp., bi-rooted) if $\left|\mathrm{Rt}\right(G\left)\right|=1$ (resp., $\left|\mathrm{Rt}\right(G\left)\right|=2$), where we call Gunrooted if $\mathrm{Rt}\left(G\right)=\varnothing$.

For a graph G, possibly with roots, a leaf-vertex is defined to be a non-root vertex $v\in V\left(G\right)\backslash \mathrm{Rt}\left(G\right)$ with degree 1, call the edge $uv$ incident to a leaf vertex v a leaf-edge, and denote $V_{\mathrm{leaf}}\left(G\right)$ and $E_{\mathrm{leaf}}\left(G\right)$ the sets of leaf-vertices and leaf-edges in G, respectively. For a graph or a rooted graph G, we define graphs $G_i,i\in {\mathbb{Z}}_{+}$ obtained from G by removing the set of leaf-vertices i times so that

$$G_0:=G;G_{i+1}:=G_i-V_{\mathrm{leaf}}\left(G_i\right),$$

where we call a vertex $v\in V_{\mathrm{leaf}}\left(G_k\right)$ a leaf k-branch and we say that a vertex $v\in V_{\mathrm{leaf}}\left(G_k\right)$ has height height ht($v)=k$ in G. The height ht($T)$ of a rooted tree T is defined to be the maximum of ht($v)$ of a vertex $v\in V\left(T\right)$. For an integer $k\ge 0$, we call a rooted tree T k-lean if T has at most one leaf k-branch. For an unrooted cyclic graph G, we regard the set of non-core-edges in G induces a collection $\mathcal{T}$ of trees each of which is rooted at a core-vertex, where we call G k-lean if each of the rooted trees in $\mathcal{T}$ is k-lean. Nearly 97% of cyclic chemical compounds with up to 100 non-hydrogen atoms in PubChem are 2-lean [24].

Two-layered Model. Let G be an unrooted graph. For an integer $\rho \ge 0$, which we call a branch-parameter, a two-layered model of G is a partition of G into an “interior” and an “exterior” in the following way. We call a vertex $v\in V\left(G\right)$ (resp., an edge $e\in E\left(G\right))$ of G an exterior-vertex (resp., exterior-edge) if ht($v)<\rho$ (resp., e is incident to an exterior-vertex) and denote the sets of exterior-vertices and exterior-edges by $V^{\mathrm{ex}}\left(G\right)$ and $E^{\mathrm{ex}}\left(G\right)$, respectively and denote $V^{\mathrm{int}}\left(G\right)=V\left(G\right)\backslash V^{\mathrm{ex}}\left(G\right)$ and $E^{\mathrm{int}}\left(G\right)=E\left(G\right)\backslash E^{\mathrm{ex}}\left(G\right)$, respectively. We call a vertex in $V^{\mathrm{int}}\left(G\right)$ (resp., an edge in $E^{\mathrm{int}}\left(G\right)$) an interior-vertex (resp., interior-edge). The set $E^{\mathrm{ex}}\left(G\right)$ of exterior-edges forms a collection of connected graphs each of which is regarded as a rooted tree T rooted at the vertex $v\in V\left(T\right)$ with the maximum ht($v)$, where we call T a ρ-fringe-tree (or a fringe-tree). Let ${\mathcal{T}}^{\mathrm{ex}}\left(G\right)$ denote the set of fringe-trees in G. The interior of G is defined to be the subgraph $(V^{\mathrm{int}}\left(G\right),E^{\mathrm{int}}\left(G\right))$ of G. Note that every core-vertex (resp., core-edge) in G is an interior-vertex (resp., interior-edge) of G. Figure 2 illustrates an example of a graph G, such that $V^{\mathrm{int}}=\{u_1,u_2,\dots ,u_{28}\}$, $V^{\mathrm{ex}}=\{w_1,w_2,\dots ,w_{19}\}$ and ${\mathcal{T}}^{\mathrm{ex}}\left(G\right)=\{T_1,T_2,\dots ,T_8\}$ for a branch-parameter $\rho =2$.

2.1.1. Modeling of Chemical Compounds

To represent a chemical compound, we assume that each chemical element $\mathrm{a}$ has a unique valence $\mathrm{val}\left(\mathrm{a}\right)\in [1,4]$ and we use a hydrogen-suppressed model, because hydrogen atoms can be added at the final stage under the assumption. In the hydrogen-suppressed model, a chemical compound C is represented by a tuple $G=(H,\alpha ,\beta )$ of a simple, connected undirected graph H and functions $\alpha :V\left(H\right)\to \mathsf{\Lambda }$ and $\beta :E\left(H\right)\to [1,3]$, where $\mathsf{\Lambda }$ is a set of non-hydrogen chemical elements such as C (carbon), O (oxygen), N (nitrogen), and so on. The set of atoms and the set of bonds in the compound C are represented by the vertex set $V\left(H\right)$ and the edge set $E\left(H\right)$, respectively. The chemical element assigned to a vertex $v\in V\left(H\right)$ is represented by $\alpha \left(v\right)$ and the bond-multiplicity between two adjacent vertices $u,v\in V\left(H\right)$ is represented by $\beta \left(e\right)$ of the edge $e=uv\in E\left(H\right)$. We say that two tuples $(H_i,{\alpha }_i,{\beta }_i),i=1,2$ are isomorphic if they admit an isomorphism $\varphi$, i.e., a bijection $\varphi :V\left(H_1\right)\to V\left(H_2\right)$ such that $uv\in E\left(H_1\right),{\alpha }_1\left(u\right)=\mathrm{a},{\alpha }_1\left(v\right)=\mathrm{b},{\beta }_1\left(uv\right)=m$ ↔ $\varphi \left(u\right)\varphi \left(v\right)\in E\left(H_2\right),{\alpha }_2\left(\varphi \left(u\right)\right)=\mathrm{a},{\alpha }_2\left(\varphi \left(v\right)\right)=\mathrm{b},{\beta }_2\left(\varphi \left(u\right)\varphi \left(v\right)\right)=m$. When $H_i$ is rooted at a vertex $r_i,i=1,2$, $(H_i,{\alpha }_i,{\beta }_i),i=1,2$ are rooted-isomorphic (r-isomorphic) if they admit an isomorphism $\varphi$ such that $\varphi \left(r_1\right)=r_2$. Chemical rooted trees $T_1$ and $T_5$ in Figure 2 are r-isomorphic.

Associated with the two functions $\alpha$ and $\beta$ in a tuple $G=(H,\alpha ,\beta )$, we introduce the following functions: ${\beta }_G:V\left(H\right)\to [0,12]$, $\mathrm{ac}:V\left(E\right)\to \mathsf{\Lambda }\times \mathsf{\Lambda }\times [1,3]$, $\mathrm{cs}:V\left(E\right)\to \mathsf{\Lambda }\times [1,4]$, and $\mathrm{ec}:V\left(E\right)\to (\mathsf{\Lambda }\times [1,4\left]\right)\times (\mathsf{\Lambda }\times [1,4\left]\right)\times [1,3]$.

For a notational convenience, we use a function ${\beta }_G:V\left(H\right)\to [0,4]$ such that ${\beta }_G\left(u\right)$ means the sum of bond-multiplicities of edges incident to a vertex u, i.e.,

$${\beta }_G\left(u\right)\triangleq \sum _{uv\in E\left(H\right)}\beta \left(uv\right)\mathrm{for}\mathrm{each}\mathrm{vertex}u\in V\left(H\right).$$

A chemical graph G is defined to be a tuple $(H,\alpha ,\beta )$ such that the valence condition at each vertex $v\in V\left(H\right)$ is satisfied, i.e.,

$${\beta }_G\left(v\right)\le \mathrm{val}\left(\alpha \left(v\right)\right),$$

where we define the hydro-degree ${deg}_{\mathrm{hyd}}\left(v\right)$ of a vertex v to be $\mathrm{val}\left(\alpha \left(v\right)\right)-{\beta }_G\left(v\right)$.

Figure 2 illustrates an example of a chemical graph $G=(H,\alpha ,\beta )$.

To represent a feature of an edge $e=uv\in E\left(H\right)$ such that $\alpha \left(u\right)=\mathrm{a}$, $\alpha \left(v\right)=\mathrm{b}$ and $\beta \left(e\right)=m$ in a chemical graph $G=(H,\alpha ,\beta )$, we use a tuple $(\mathrm{a},\mathrm{b},m)\in \mathsf{\Lambda }\times \mathsf{\Lambda }\times [1,3]$, which we call the adjacency-configuration $\mathrm{ac}\left(e\right)$ of the edge e. We introduce a total order < over the elements in $\mathsf{\Lambda }$ to distinguish with $(\mathrm{a},\mathrm{b},m)$ and $(\mathrm{b},\mathrm{a},m)$ $(\mathrm{a}\ne \mathrm{b})$ notationally. For a tuple $\nu =(\mathrm{a},\mathrm{b},m)$, let $\overline{\nu }$ denote the tuple $(\mathrm{b},\mathrm{a},m)$.

To represent a feature of a vertex $v\in V\left(H\right)$ with $\alpha \left(v\right)=\mathrm{a}$ that has d atoms in its neighbor in a chemical graph $G=(H,\alpha ,\beta )$, we use a pair $(\mathrm{a},d)\in \mathsf{\Lambda }\times [1,4]$, which we call the chemical symbol $\mathrm{cs}\left(v\right)$ of the vertex v. We treat $(\mathrm{a},d)$ as a single symbol $\mathrm{a}d$, and define ${\mathsf{\Lambda }}_{\mathrm{dg}}$ to be the set of all chemical symbols $\mu =\mathrm{a}d\in \mathsf{\Lambda }\times [1,4]$.

To represent a feature of an edge $e=uv\in E\left(H\right)$ such that $\mathrm{cs}\left(u\right)=\mu$, $\mathrm{cs}\left(v\right)=\xi$ and $\beta \left(e\right)=m$ in a chemical graph $G=(H,\alpha ,\beta )$, we use a tuple $(\mu ,\xi ,m)\in {\mathsf{\Lambda }}_{\mathrm{dg}}\times {\mathsf{\Lambda }}_{\mathrm{dg}}\times [1,3]$, which we call the edge-configuration $\mathrm{ec}\left(e\right)$ of the edge e. We introduce a total order < over the elements in ${\mathsf{\Lambda }}_{\mathrm{dg}}$ to distinguish with $(\mu ,\xi ,m)$ and $(\xi ,\mu ,m)$ $(\mu \ne \xi )$ notationally. For a tuple $\gamma =(\mu ,\xi ,m)$, let $\overline{\gamma }$ denote the tuple $(\xi ,\mu ,m)$.

To represent a feature of the exterior of a chemical graph $G=(H,\alpha ,\beta )$, a $\rho$-fringe-tree in ${\mathcal{T}}^{\mathrm{ex}}\left(G\right)$ is called a fringe-configuration in the exterior.

2.1.2. Introducing Descriptors of Feature Vectors

This section introduces descriptors to define our feature vectors. Let $\pi$ be a chemical property for which we will construct a prediction function ${\eta }_{\mathcal{N}}$ from a feature vector $f\left(G\right)$ of a chemical graph to a predicted value $y\in \mathbb{R}$ for the chemical property of G.

We first choose a set $\mathsf{\Lambda }$ of non-hydrogen chemical elements and then collect a data set $D_{\pi }$ of chemical compounds C whose chemical elements belong to $\mathsf{\Lambda }\cup \left\{\mathrm{H}\right\}$, where we regard $D_{\pi }$ as a set of chemical graphs that represent the chemical compounds C in $D_{\pi }$. To define the interior/exterior of chemical graphs $G\in D_{\pi }$, we next choose a branch-parameter $\rho$, where we recommend $\rho =2$.

Let ${\mathsf{\Lambda }}^{\mathrm{int}}\left(D_{\pi }\right)$ (resp., ${\mathsf{\Lambda }}^{\mathrm{ex}}\left(D_{\pi }\right)$) denote the set of chemical elements used in the set of interior-vertices (resp., exterior-vertices) over all chemical graphs $G\in D_{\pi }$, and ${\mathsf{\Gamma }}^{\mathrm{int}}\left(D_{\pi }\right)$ denote the set of edge-configurations used in the set of interior-edges over all chemical graphs $G\in D_{\pi }$. Let $\mathcal{F}\left(D_{\pi }\right)$ denote the set of chemical rooted trees $\psi$ r-isomorphic to a $\rho$-fringe-tree $T\in {\mathcal{T}}^{\mathrm{ex}}\left(G\right)$ over all chemical graphs $G\in D_{\pi }$.

We define an integer encoding of a finite set A of elements to be a bijection $\sigma :A\to [1,|A\left|\right]$, where we denote by $\left[A\right]$ the set $[1,|A\left|\right]$ of integers. Introduce an integer coding of each of the sets ${\mathsf{\Lambda }}^{\mathrm{int}}\left(D_{\pi }\right)$, ${\mathsf{\Lambda }}^{\mathrm{ex}}\left(D_{\pi }\right)$, ${\mathsf{\Gamma }}^{\mathrm{int}}\left(D_{\pi }\right)$ and $\mathcal{F}\left(D_{\pi }\right)$. Let ${\left[\mathrm{a}\right]}^{\mathrm{int}}$ (resp., ${\left[\mathrm{a}\right]}^{\mathrm{ex}}$) denote the coded integer of an element $\mathrm{a}\in {\mathsf{\Lambda }}^{\mathrm{int}}\left(D_{\pi }\right)$ (resp., $\mathrm{a}\in {\mathsf{\Lambda }}^{\mathrm{ex}}\left(D_{\pi }\right)$), $\left[\gamma \right]$ denote the coded integer of an element $\gamma$ in ${\mathsf{\Gamma }}^{\mathrm{int}}\left(D_{\pi }\right)$ and $\left[\psi \right]$ denote an element $\psi$ in $\mathcal{F}\left(D_{\pi }\right)$.

For each chemical element $\mathrm{a}\in \mathsf{\Lambda }$, let $\mathrm{mass}\left(\mathrm{a}\right)$ and $\mathrm{val}\left(\mathrm{a}\right)$ denote the mass and valence of $\mathrm{a}$, respectively. In our model, we use integers ${\mathrm{mass}}^{*}\left(\mathrm{a}\right)=\lfloor 10·\mathrm{mass}\left(\mathrm{a}\right)\rfloor$, $\mathrm{a}\in \mathsf{\Lambda }$.

We define the feature vector $f\left(G\right)$ of a chemical graph $G=(H,\alpha ,\beta )\in D_{\pi }$ to be a vector that consists of the following non-negative integer descriptors ${\mathrm{dcp}}_i\left(G\right)$, $i\in [1,K]$, where $K=17+|{\mathsf{\Lambda }}^{\mathrm{int}}\left(D_{\pi }\right)|+|{\mathsf{\Lambda }}^{\mathrm{ex}}\left(D_{\pi }\right)|+|{\mathsf{\Gamma }}^{\mathrm{int}}\left(D_{\pi }\right)|+|\mathcal{F}\left(D_{\pi }\right)|$.

• ${\mathrm{dcp}}_1\left(G\right)$: the number $n\left(G\right)=\left|V\right(G\left)\right|$ of vertices in G.
• ${\mathrm{dcp}}_2\left(G\right)$: the number $|V^{\mathrm{int}}\left(G\right)|$ of interior-vertices in G.
• ${\mathrm{dcp}}_3\left(G\right)$: the average $\overline{\mathrm{ms}}\left(G\right)$ of mass${}^{*}$ over all non-hydrogen atoms in G, i.e., $\overline{\mathrm{ms}}\left(G\right)\triangleq {\sum }_{v\in V\left(G\right)}{\mathrm{mass}}^{*}\left(\alpha \left(v\right)\right)/n\left(G\right)$.
• ${\mathrm{dcp}}_i\left(G\right)$, $i=3+d,d\in [1,4]$: the number ${\mathrm{dg}}_d\left(G\right)$ of interior-vertices of degree d in G.
• ${\mathrm{dcp}}_i\left(G\right)$, $i=7+d,d\in [1,4]$: the number ${\mathrm{dg}}_d^{\mathrm{int}}\left(G\right)$ of interior-vertices of interior-degree ${deg}_{(V^{\mathrm{int}},E^{\mathrm{int}})}\left(v\right)=d$ in the interior $(V^{\mathrm{int}},E^{\mathrm{int}})$ of G.
• ${\mathrm{dcp}}_i\left(G\right)$, $i=11+d,d\in [0,3]$: the number ${\mathrm{hydg}}_d\left(G\right)$ of vertices in G of hydro-degree ${deg}_{\mathrm{hyd}}\left(v\right)=d$.
• ${\mathrm{dcp}}_i\left(G\right)$, $i=15+m$, $m\in [2,3]$: the number ${\mathrm{bd}}_m^{\mathrm{int}}\left(G\right)$ of interior-edges with bond multiplicity m in G, i.e., ${\mathrm{bd}}_m^{\mathrm{int}}\left(G\right)\triangleq \{e\in E^{\mathrm{int}}\mid \beta \left(e\right)=m\}$.
• ${\mathrm{dcp}}_i\left(G\right)$, $i=17+{\left[\mathrm{a}\right]}^{\mathrm{int}}$, $\mathrm{a}\in {\mathsf{\Lambda }}^{\mathrm{int}}\left(D_{\pi }\right)$: the frequency ${\mathrm{na}}_{\mathrm{a}}^{\mathrm{int}}\left(G\right)$ of chemical element $\mathrm{a}$ in the set of interior-vertices in G.
• ${\mathrm{dcp}}_i\left(G\right)$, $i=17+|{\mathsf{\Lambda }}^{\mathrm{int}}\left(D_{\pi }\right)|+{\left[\mathrm{a}\right]}^{\mathrm{ex}}$, $\mathrm{a}\in {\mathsf{\Lambda }}^{\mathrm{ex}}\left(D_{\pi }\right)$: the frequency ${\mathrm{na}}_{\mathrm{a}}^{\mathrm{ex}}\left(G\right)$ of chemical element $\mathrm{a}$ in the set of exterior-vertices in G.
• ${\mathrm{dcp}}_i\left(G\right)$, $i=17+|{\mathsf{\Lambda }}^{\mathrm{int}}\left(D_{\pi }\right)|+|{\mathsf{\Lambda }}^{\mathrm{ex}}\left(D_{\pi }\right)|+\left[\gamma \right]$, $\gamma \in {\mathsf{\Gamma }}^{\mathrm{int}}\left(D_{\pi }\right)$: the frequency ${\mathrm{ec}}_{\gamma }\left(G\right)$ of edge-configuration $\gamma$ in the set of interior-edges $e\in E^{\mathrm{int}}$ in G.
• ${\mathrm{dcp}}_i\left(G\right)$, $i=17+|{\mathsf{\Lambda }}^{\mathrm{int}}\left(D_{\pi }\right)|+|{\mathsf{\Lambda }}^{\mathrm{ex}}\left(D_{\pi }\right)|+|{\mathsf{\Gamma }}^{\mathrm{int}}\left(D_{\pi }\right)|+\left[\psi \right]$, $\psi \in \mathcal{F}\left(D_{\pi }\right)$: the frequency ${\mathrm{fc}}_{\psi }\left(G\right)$ of fringe-configuration $\psi$ in the set of $\rho$-fringe-trees in G.

2.2. Specifying Target Chemical Graphs

Given a prediction function ${\eta }_{\mathcal{N}}$ and a target value $y^{*}\in \mathbb{R}$, we call a chemical graph $G^{*}$ such that ${\eta }_{\mathcal{N}}\left(x^{*}\right)=y^{*}$ for the feature vector $x^{*}=f\left(G^{*}\right)$ a target chemical graph. This section presents a set of rules for specifying topological substructure of a target chemical graph in a flexible way in Stage 4.

We first describe how to reduce a chemical graph $G=(H,\alpha ,\beta )$ into an abstract form based on which our specification rules will be defined. To illustrate the reduction process, we use the chemical graph $G=(H,\alpha ,\beta )$ in Figure 2.

R1

Removal of all $\rho$-fringe-trees: The interior $H^{\mathrm{int}}=(V^{\mathrm{int}}\left(H\right),E^{\mathrm{int}}\left(H\right))$ of G is obtained by removing the non-root vertices of each $\rho$-fringe-trees $T\in {\mathcal{T}}^{\mathrm{ex}}\left(G\right)$. Figure 4 illustrates the interior $H^{\mathrm{int}}$ of chemical graph G with $\rho =2$ in Figure 2.

R2

Removal of some leaf paths: We call a $u,v$-path Q in $H^{\mathrm{int}}$ a leaf path if vertex v is a leaf-vertex of $H^{\mathrm{int}}$ and the degree of each internal vertex of Q in $H^{\mathrm{int}}$ is 2, where we regard that Q is rooted at vertex u. A connected subgraph S of the interior $H^{\mathrm{int}}$ of G is called a cyclical-base if S is obtained from H by removing the vertices in $V\left(Q_u\right)\backslash \left\{u\right\},u\in X$ for a subset X of interior-vertices and a set $\{Q_u\mid u\in X\}$ of leaf $u,v$-paths $Q_u$ such that no two paths $Q_u$ and $Q_{u^{\prime }}$ share a vertex. Figure 5a illustrates a cyclical-base $S=H^{\mathrm{int}}-{\bigcup }_{u\in X}(V\left(Q_u\right)\backslash \left\{u\right\})$ of the interior $H^{\mathrm{int}}$ for a set $\{Q_{u_5}=(u_5,u_{24}),Q_{u_{18}}=(u_{18},u_{25},u_{26},u_{27}),Q_{u_{22}}=(u_{22},u_{28})\}$ of leaf paths in Figure 4.

R3

Contraction of some pure paths: A path in S is called pure if each internal vertex of the path is of degree 2. Choose a set $\mathcal{P}$ of several pure paths in S so that no two paths share vertices except for their end-vertices. A graph $S^{\prime }$ is called a contraction of a graph S (with respect to $\mathcal{P}$) if $S^{\prime }$ is obtained from S by replacing each pure $u,v$-path with a single edge $a=uv$, where $S^{\prime }$ may contain multiple edges between the same pair of adjacent vertices. Figure 5b illustrates a contraction $S^{\prime }$ obtained from the chemical graph S by contracting each $uv$-path $P_a\in \mathcal{P}$ into a new edge $a=uv$, where $a_1=u_1{u}_2,a_2=u_1{u}_3,a_3=u_4{u}_7,a_4=u_{10}{u}_{11}$, and $a_5=u_{11}{u}_{12}$, and $\mathcal{P}=\{P_{a_1}=(u_1,u_{13},u_2),P_{a_2}=(u_1,u_{14},u_3),P_{a_3}=(u_4,u_{15},u_{16},u_7),P_{a_4}=(u_{10},u_{17},u_{18},u_{19},u_{11}),P_{a_5}=(u_{11},u_{20},u_{21},u_{22},u_{12})\}$ of pure paths in Figure 5a.

We will define a set of rules so that a chemical graph can be obtained from a graph (called a seed graph in the next section) by applying processes R3 to R1 in a reverse way. We specify topological substructures of a target chemical graph with a tuple $(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ called a target specification defined under the set of the following rules.

Seed Graphs

A seed graph $G_{\mathrm{C}}=(V_{\mathrm{C}},E_{\mathrm{C}})$ is defined to be a graph (possibly with multiple edges) such that the edge set $E_{\mathrm{C}}$ consists of four sets $E_{(\ge 2)}$, $E_{(\ge 1)}$, $E_{(0/1)}$, and $E_{(=1)}$, where each of them can be empty. A seed graph plays a role of the most abstract form $S^{\prime }$ in R3. Figure 3a illustrates an example of a seed graph, where $V_{\mathrm{C}}=\{u_1,u_2,\dots ,u_{12}\}$, $E_{(\ge 2)}=\{a_1,a_2,\dots ,a_5\}$, $E_{(\ge 1)}=\left\{a_6\right\}$, $E_{(0/1)}=\left\{a_7\right\}$, and $E_{(=1)}=\{a_8,a_9,\dots ,a_{17}\}$.

A subdivision S of $G_{\mathrm{C}}$ is a graph constructed from a seed graph $G_{\mathrm{C}}$ according to the following rules:

-

Each edge $e=uv\in E_{(\ge 2)}$ is replaced with a $u,v$-path $P_e$ of length at least 2;

-

Each edge $e=uv\in E_{(\ge 1)}$ is replaced with a $u,v$-path $P_e$ of length at least 1 (equivalently e is directly used or replaced with a $u,v$-path $P_e$ of length at least 2);

-

Each edge $e\in E_{(0/1)}$ is either used or discarded; and

-

Each edge $e\in E_{(=1)}$ is always used directly.

We allow a possible elimination of edges in $E_{(0/1)}$ as an optional rule in constructing a target chemical graph from a seed graph, even though such an operation has not been included in the process R3. A subdivision S plays a role of a cyclical-base in R2. A target chemical graph $G=(H,\alpha ,\beta )$ will contain S as a subgraph of the interior $H^{\mathrm{int}}$ of G.

Interior-Specification

A graph $H^{*}$ that serves as the interior $H^{\mathrm{int}}$ of a target chemical graph G will be constructed as follows. First, construct a subdivision S of a seed graph $G_{\mathrm{C}}$ by replacing each edge edge $e=u{u}^{\prime }\in E_{(\ge 2)}\cup E_{(\ge 1)}$ with a pure $u,u^{\prime }$-path $P_e$. Next, construct a supergraph $H^{*}$ of S by attaching a leaf path $Q_v$ at each vertex $v\in V_{\mathrm{C}}$ or at an internal vertex $v\in V\left(P_e\right)\backslash \{u,u^{\prime }\}$ of each pure $u,u^{\prime }$-path $P_e$ for some edge $e=u{u}^{\prime }\in E_{(\ge 2)}\cup E_{(\ge 1)}$, where possibly $Q_v=v,E\left(Q_v\right)=\varnothing$ (i.e., we do not attach any new edges to v). We introduce the following rules for specifying the size of $H^{*}$, the length $\left|E\right(P_e\left)\right|$ of a pure path $P_e$, the length $\left|E\right(Q_v\left)\right|$ of a leaf path $Q_v$, the number of leaf paths $Q_v$, and a bond-multiplicity of each interior-edge, where we call the set of prescribed constants an interior-specification ${\sigma }_{\mathrm{int}}$:

-

Lower and upper bounds ${\mathrm{n}}_{\mathrm{LB}}^{\mathrm{int}},{\mathrm{n}}_{\mathrm{UB}}^{\mathrm{int}}\in {\mathbb{Z}}_{+}$ on the number of interior-vertices of a target chemical graph G.

-

For each edge $e=u{u}^{\prime }\in E_{(\ge 2)}\cup E_{(\ge 1)}$,

• a lower bound ${\ell }_{\mathrm{LB}}\left(e\right)$ and an upper bound ${\ell }_{\mathrm{UB}}\left(e\right)$ on the length $\left|E\right(P_e\left)\right|$ of a pure $u,u^{\prime }$-path $P_e$. (For a notational convenience, set ${\ell }_{\mathrm{LB}}\left(e\right):=0$, ${\ell }_{\mathrm{UB}}\left(e\right):=1$, $e\in E_{(0/1)}$ and ${\ell }_{\mathrm{LB}}\left(e\right):=1$, ${\ell }_{\mathrm{UB}}\left(e\right):=1$, $e\in E_{(=1)}$. )
• a lower bound ${\mathrm{bl}}_{\mathrm{LB}}\left(e\right)$ and an upper bound ${\mathrm{bl}}_{\mathrm{UB}}\left(e\right)$ on the number of leaf paths $Q_v$ attached to at internal vertices v of a pure $u,u^{\prime }$-path $P_e$.
• a lower bound ${\mathrm{ch}}_{\mathrm{LB}}\left(e\right)$ and an upper bound ${\mathrm{ch}}_{\mathrm{UB}}\left(e\right)$ on the maximum length $\left|E\right(Q_v\left)\right|$ of a leaf path $Q_v$ attached at an internal vertex $v\in V\left(P_e\right)\backslash \{u,u^{\prime }\}$ of a pure $u,u^{\prime }$-path $P_e$.

-

For each vertex $v\in V_{\mathrm{C}}$,

• a lower bound ${\mathrm{ch}}_{\mathrm{LB}}\left(e\right)$ and an upper bound ${\mathrm{ch}}_{\mathrm{UB}}\left(e\right)$ on the number of leaf paths $Q_v$ attached to v, where $0\le {\mathrm{ch}}_{\mathrm{LB}}\left(e\right)\le {\mathrm{ch}}_{\mathrm{UB}}\left(e\right)\le 1$.
• a lower bound ${\mathrm{ch}}_{\mathrm{LB}}\left(v\right)$ and an upper bound ${\mathrm{ch}}_{\mathrm{UB}}\left(v\right)$ on the length $\left|E\right(Q_v\left)\right|$ of a leaf path $Q_v$ attached to v.

-

For each edge $e=u{u}^{\prime }\in E_{\mathrm{C}}$, a lower bound ${\mathrm{bd}}_{m,\mathrm{LB}}\left(e\right)$ and an upper bound ${\mathrm{bd}}_{m,\mathrm{UB}}\left(e\right)$ on the number of edges with bond-multiplicity $m\in [2,3]$ in $u,u^{\prime }$-path $P_e$, where we regard $P_e$, $e\in E_{(0/1)}\cup E_{(=1)}$ as single edge e.

We call a graph $H^{*}$ that satisfies an interior-specification ${\sigma }_{\mathrm{int}}$ a ${\sigma }_{\mathrm{int}}$-extension of $G_{\mathrm{C}}$, where the bond-multiplicity of each edge has been determined.

Table 1. Example 1 of an interior-specification ${\sigma }_{\mathrm{int}}$.

shows an example of an interior-specification ${\sigma }_{\mathrm{int}}$ to the seed graph $G_{\mathrm{C}}$ in Figure 3.

Figure 6 illustrates an example of an ${\sigma }_{\mathrm{int}}$-extension $H^{*}$ of seed graph $G_{\mathrm{C}}$ in Figure 3 under the interior-specification ${\sigma }_{\mathrm{int}}$ in Table 1.

Chemical-Specification

Let $H^{*}$ be a graph that serves as the interior $H^{\mathrm{int}}$ of a target chemical graph G, where the bond-multiplicity of each edge in $H^{*}$ has be determined. Finally, we introduce a set of rules for constructing a target chemical graph G from $H^{*}$ by choosing a chemical element $\mathrm{a}\in \mathsf{\Lambda }$ and assigning a $\rho$-fringe-tree $\psi$ to each interior-vertex $v\in V^{\mathrm{int}}$. We introduce the following rules for specifying the size of G, a set of chemical rooted trees that are allowed to use as $\rho$-fringe-trees and lower and upper bounds on the frequency of a chemical element, a chemical symbol, and an edge-configuration, where we call the set of prescribed constants a chemical specification ${\sigma }_{\mathrm{ce}}$:

-

Lower and upper bounds $n_{\mathrm{LB}},n^{*}\in {\mathbb{Z}}_{+}$ on the number of vertices in G, where ${\mathrm{n}}_{\mathrm{LB}}^{\mathrm{int}}\le n_{\mathrm{LB}}\le n^{*}$.

-

Subsets $\mathcal{F}\left(v\right)\subseteq \mathcal{F}\left(D_{\pi }\right),v\in V_{\mathrm{C}}$ and ${\mathcal{F}}_E\subseteq \mathcal{F}\left(D_{\pi }\right)$ of chemical rooted trees with height at most $\rho$, where we require that every $\rho$-fringe-tree $T_v$ rooted at a vertex $v\in V_{\mathrm{C}}$ (resp., at an internal vertex v not in $V_{\mathrm{C}}$) in G belongs to $\mathcal{F}\left(v\right)$ (resp., ${\mathcal{F}}_E$). Let ${\mathcal{F}}^{*}:={\mathcal{F}}_E\cup {\bigcup }_{v\in V_{\mathrm{C}}}\mathcal{F}\left(v\right)$ and ${\mathsf{\Lambda }}^{\mathrm{ex}}$ denote the set of chemical elements assigned to non-root vertices over all chemical rooted trees in ${\mathcal{F}}^{*}$.

-

A subset ${\mathsf{\Lambda }}^{\mathrm{int}}\subseteq {\mathsf{\Lambda }}^{\mathrm{int}}\left(D_{\pi }\right)$, where we require that every chemical element $\alpha \left(v\right)$ assigned to an interior-vertex v in G belongs to ${\mathsf{\Lambda }}^{\mathrm{int}}$. Let $\mathsf{\Lambda }:={\mathsf{\Lambda }}^{\mathrm{int}}\cup {\mathsf{\Lambda }}^{\mathrm{ex}}$ and ${\mathrm{na}}_{\mathrm{a}}\left(G\right)$ (resp., ${\mathrm{na}}_{\mathrm{a}}^{\mathrm{int}}\left(G\right)$ and ${\mathrm{na}}_{\mathrm{a}}^{\mathrm{ex}}\left(G\right)$) denote the number of vertices (resp., interior-vertices and exterior-vertices) v such that $\alpha \left(v\right)=\mathrm{a}$ in G.

-

A set ${\mathsf{\Lambda }}_{\mathrm{dg}}^{\mathrm{int}}\subseteq \mathsf{\Lambda }\times [1,4]$ of chemical symbols and a set ${\mathsf{\Gamma }}^{\mathrm{int}}\subseteq {\mathsf{\Gamma }}^{\mathrm{int}}\left(D_{\pi }\right)$ of edge-configurations $(\mu ,\xi ,m)$ with $\mu \le \xi$, where we require that the edge-configuration $\mathrm{ec}\left(e\right)$ of an interior-edge e in G belongs to ${\mathsf{\Gamma }}^{\mathrm{int}}$. We do not distinguish $(\mu ,\xi ,m)$ and $(\xi ,\mu ,m)$.

-

Define ${\mathsf{\Gamma }}_{\mathrm{ac}}^{\mathrm{int}}$ to be the set of adjacency-configurations such that ${\mathsf{\Gamma }}_{\mathrm{ac}}^{\mathrm{int}}:=\{(\mathrm{a},\mathrm{b},m)\mid (\mathrm{a}d,\mathrm{b}{d}^{\prime },m)\in {\mathsf{\Gamma }}^{\mathrm{int}}\}$. Let ${\mathrm{ac}}_{\nu }^{\mathrm{int}}\left(G\right),\nu \in {\mathsf{\Gamma }}_{\mathrm{ac}}^{\mathrm{int}}$ denote the number of interior-edges e such that $\mathrm{ac}\left(e\right)=\nu$ in G.

-

Subsets ${\mathsf{\Lambda }}^{*}\left(v\right)\subseteq \{\mathrm{a}\in {\mathsf{\Lambda }}^{\mathrm{int}}\mid \mathrm{val}\left(\mathrm{a}\right)\ge 2\}$, $v\in V_{\mathrm{C}}$, we require that every chemical element $\alpha \left(v\right)$ assigned to a vertex $v\in V_{\mathrm{C}}$ in the seed graph belongs to ${\mathsf{\Lambda }}^{*}\left(v\right)$.

-

Lower and upper bound functions ${\mathrm{na}}_{\mathrm{LB}},{\mathrm{na}}_{\mathrm{UB}}:\mathsf{\Lambda }\to [1,n^{*}]$ and ${\mathrm{na}}_{\mathrm{LB}}^{\mathrm{int}},{\mathrm{na}}_{\mathrm{UB}}^{\mathrm{int}}:{\mathsf{\Lambda }}^{\mathrm{t}}\to [1,n^{*}]$ on the number of interior-vertices v such that $\alpha \left(v\right)=\mathrm{a}$ in G.

-

Lower and upper bound functions ${\mathrm{ns}}_{\mathrm{LB}}^{\mathrm{int}},{\mathrm{ns}}_{\mathrm{UB}}^{\mathrm{int}}:{\mathsf{\Lambda }}_{\mathrm{dg}}^{\mathrm{int}}\to [1,n^{*}]$ on the number of interior-vertices v such that $\mathrm{cs}\left(v\right)=\mu$ in G.

-

Lower and upper bound functions ${\mathrm{ac}}_{\mathrm{LB}}^{\mathrm{int}},{\mathrm{ac}}_{\mathrm{UB}}^{\mathrm{int}}:{\mathsf{\Gamma }}_{\mathrm{ac}}^{\mathrm{int}}\to {\mathbb{Z}}_{+}$ on the number of interior-edges e such that $\mathrm{ac}\left(e\right)=\nu$ in G.

-

Lower and upper bound functions ${\mathrm{ec}}_{\mathrm{LB}}^{\mathrm{int}},{\mathrm{ec}}_{\mathrm{UB}}^{\mathrm{int}}:{\mathsf{\Gamma }}^{\mathrm{int}}\to {\mathbb{Z}}_{+}$ on the number of interior-edges e such that $\mathrm{ec}\left(e\right)=\gamma$ in G.

We call a chemical graph G that satisfies a chemical specification ${\sigma }_{\mathrm{ce}}$ a $({\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$-extension of $G_{\mathrm{C}}$, and denote by $\mathcal{G}(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ the set of all $({\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$-extensions of $G_{\mathrm{C}}$.

Table 2. Example 2 of a chemical-specification ${\sigma }_{\mathrm{ce}}$.

shows an example of a chemical-specification ${\sigma }_{\mathrm{ce}}$ to the seed graph $G_{\mathrm{C}}$ in Figure 3.

Figure 2 illustrates an example of a $({\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$-extension of $G_{\mathrm{C}}$ obtained from the ${\sigma }_{\mathrm{int}}$-extension $H^{*}$ in Figure 6 under the chemical-specification ${\sigma }_{\mathrm{ce}}$ in Table 2.

Our specification of topological substructures is similar to that proposed by Akutsu and Nagamochi [27], wherein a target chemical graph is restricted to $\rho$-lean cyclic graphs and prescribed substructures cannot be specified in the acyclic part. In our new method, a chemical graph with any structure can be handled and substructures in the acyclic part can be fixed.

2.3. Examples of Specification

We here present some cases where a target specification $(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ can be chosen based on a set ${\mathcal{G}}^{*}$ of given chemical graphs with a similar structure so that ${\mathcal{G}}^{*}$ becomes a subset of $\mathcal{G}(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$. In such a case, every target chemical graph in $\mathcal{G}(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ possesses a common structure over the given set ${\mathcal{G}}^{*}$.

Figure 7 illustrates a set ${\mathcal{G}}^{*}$ of four flavonoids and a seed graph $G_{\mathrm{C}}$ for $\rho =2$ so that ${\mathcal{G}}^{*}\subseteq \mathcal{G}(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ for a choice of an interior-specification ${\sigma }_{\mathrm{int}}$ and a chemical-specification ${\sigma }_{\mathrm{ce}}$. Let $\mathsf{\Lambda }:=\left\{\mathrm{C},\mathrm{O}\right\}$. In the seed graph $G_{\mathrm{C}}=(V_{\mathrm{C}},E_{\mathrm{C}})$, we set $E_{(\ge 1)}:=\{a_1,a_2\}$, $E_{(0/1)}:=\left\{a_3\right\}$, and $E_{(=1)}:=E_{\mathrm{C}}\backslash \{a_1,a_2,a_3\}$ and predetermine the chemical element $\alpha \left(u\right)$ for each vertex $u\in V_{\mathrm{C}}$ and the bond-multiplicity $\beta \left(e\right)$ for each edge $e\in E_{(=1)}$ as in Figure 7e, i.e., ${\mathsf{\Lambda }}^{*}\left(u\right):=\left\{\mathrm{a}\right\}$ for $\mathrm{a}=\alpha \left(u\right)$ and ${\mathrm{bd}}_{m,\mathrm{LB}}\left(e\right):=1$ for $m=\beta \left(e\right)$. Figure 7f illustrates a set ${\mathcal{F}}^{*}$ of chemical rooted trees for the 2-fringe-trees in a target chemical graph. For vertices in $G_{\mathrm{C}}$, we set ${\mathrm{ch}}_{\mathrm{UB}}\left(u\right):=0,u\in V_{\mathrm{C}}$, $\mathcal{F}\left(u_i\right):=\left\{{\psi }_3\right\},i\in [1,3]$, $\mathcal{F}\left(u_4\right):=\{{\psi }_1,{\psi }_3\}$, $\mathcal{F}\left(u_5\right):=\left\{{\psi }_4\right\}$, $\mathcal{F}\left(u_6\right):=\left\{{\psi }_2\right\}$, and $\mathcal{F}\left(u\right):=\left\{{\psi }_1\right\},u\in V_{\mathrm{C}}\backslash \{u_1,u_2,\dots ,u_6\}$. For edges $a_i\in E_{(\ge 1)},i=1,2$, we set ${\ell }_{\mathrm{UB}}\left(a_i\right):=2,{\mathrm{ch}}_{\mathrm{UB}}\left(a_i\right):=0$ and ${\mathcal{F}}_E:=\{{\psi }_1,{\psi }_2\}$, where a pure path $P_{a_i}$ may be introduced in a target chemical graph. We see that every given chemical graph $G_i\in {\mathcal{G}}^{*}$ belongs to $\mathcal{G}(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ by setting the other specification in ${\sigma }_{\mathrm{int}}$ and ${\sigma }_{\mathrm{ce}}$ adequately.

Figure 8 illustrates a set ${\mathcal{G}}^{*}$ of three dibenzodiazepine atypical antipsychotics, and a seed graph $G_{\mathrm{C}}$ for $\rho =2$ so that ${\mathcal{G}}^{*}\subseteq \mathcal{G}(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ for a choice of an interior-specification ${\sigma }_{\mathrm{int}}$ and a chemical-specification ${\sigma }_{\mathrm{ce}}$. Let $\mathsf{\Lambda }:=\left\{\mathrm{C},\mathrm{O},\mathrm{N},\mathrm{S},\mathrm{Cl}\right\}$. In the seed graph $G_{\mathrm{C}}=(V_{\mathrm{C}},E_{\mathrm{C}})$, we set $E_{(\ge 2)}:=\left\{a_1\right\}$ and $E_{(=1)}:=E_{\mathrm{C}}\backslash \left\{a_1\right\}$ and predetermine the chemical element $\alpha \left(u\right)$ for each vertex $u\in V_{\mathrm{C}}$ and the bond-multiplicity $\beta \left(e\right)$ for each edge $e\in E_{(=1)}$ as in Figure 8d. Figure 8e illustrates a set ${\mathcal{F}}^{*}$ of chemical rooted trees for the 2-fringe-trees in a target chemical graph. For vertices in $G_{\mathrm{C}}$, we set ${\mathrm{ch}}_{\mathrm{UB}}\left(u\right):=0,u\in V_{\mathrm{C}}\backslash \left\{u_2\right\}$, ${\mathrm{ch}}_{\mathrm{LB}}\left(u_2\right):=0$, ${\mathrm{ch}}_{\mathrm{UB}}\left(u_2\right):=4$, $\mathcal{F}\left(u_1\right):=\{{\psi }_3,{\psi }_7\}$, $\mathcal{F}\left(u_2\right):=\{{\psi }_1,{\psi }_6\}$, $\mathcal{F}\left(u_i\right):=\left\{{\psi }_3\right\},i\in [3,5]$, and $\mathcal{F}\left(u\right):=\left\{{\psi }_1\right\},u\in V_{\mathrm{C}}\backslash \{u_1,u_2,\dots ,u_6\}$, where a leaf path $Q_{u_2}$ may be introduced in a target chemical graph. For edge $a_1\in E_{(\ge 2)}$, we set ${\ell }_{\mathrm{UB}}\left(a_1\right):=3,{\mathrm{ch}}_{\mathrm{UB}}\left(a_1\right):=0$ and ${\mathcal{F}}_E:=\{{\psi }_1,{\psi }_2,{\psi }_4,{\psi }_8\}$. We see that every given chemical graph $G_i\in {\mathcal{G}}^{*}$ belongs to $\mathcal{G}(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ by setting the other specification in ${\sigma }_{\mathrm{int}}$ and ${\sigma }_{\mathrm{ce}}$ adequately.

3. Results

We implemented our method of Stages 1 to 5 for inferring chemical graphs under a given target specification and conducted experiments to evaluate the computational efficiency. We executed the experiments on a PC with Processor: 3.0 GHz Core i7-9700 (3.0 GHz) Memory: 16 GB RAM DDR4. We used ChemDoodle version 10.2.0 for constructing 2D drawings of chemical graphs.

To conduct experiments for Stages 1 to 5, we selected six chemical properties $\pi$: octanol/water partition coefficient (K_OW), boiling point (B_P), melting point (M_P), flash point (closed cup) (F_P), lipophylicity (L_P), solubility (S_L) provided by HSDB from PubChem [29] for K_OW, B_P, M_P, and F_P, figshare [30] for L_P and MoleculeNet [31] for S_L.

Results on Phase 1.

We implemented Stages 1, 2, and 3 in Phase 1 as follows.

Stage 1. We set a graph class $\mathcal{G}$ to be the set of all chemical graphs with any graph structure, and set a branch-parameter $\rho$ to be 2. For each property $\pi \in \{$K_OW, B_P, M_P, F_P, L_P, S_L}, we first select a set $\mathsf{\Lambda }$ of chemical elements and then collect a data set $D_{\pi }$ on chemical graphs over the set $\mathsf{\Lambda }$ of chemical elements. To construct the data set $D_{\pi }$, we eliminated chemical compounds that have at most three carbon atoms or contain a charged element such as ${\mathrm{N}}^{+}$ or an element $\mathrm{a}\in \mathsf{\Lambda }$ whose valence is different from our setting of valence function $\mathrm{val}$.

Table 3. Data sets for stage 1 in phase 1.

$\mathit{\pi }$	$\mathsf{\Lambda }$	$|{\mathit{D}}_{\mathit{\pi }}|$	$|{\mathsf{\Gamma }}^{\mathbf{int}}\left({\mathit{D}}_{\mathit{\pi }}\right)|$	$\left|\mathcal{F}\right({\mathit{D}}_{\mathit{\pi }}\left)\right|$	$[\underline{\mathit{n}},\overline{\mathit{n}}]$	$[\underline{\mathit{a}},\overline{\mathit{a}}]$
KOW	C,O,N	644	24	109	[4, 58]	[−7.53, 13.45]
KOW	C,O,N,S,Cl	837	31	142	[4, 73]	[−7.53, 13.45]
BP	C,O,N	358	21	91	[4, 30]	[−11.70, 470.0]
BP	C,O,N,S,Cl	425	23	114	[4, 30]	[−11.70, 470.0]
MP	C,O,N	448	22	94	[4, 122]	[−185.3, 300.0]
MP	C,O,N,S,Cl	548	26	118	[4, 122]	[−185.3, 300.0]
FP	C,O,N	348	20	85	[4, 66]	[−82.99, 300.0]
FP	C,O,N,S,Cl	399	24	107	[4, 66]	[−82.99, 300.0]
LP	C,O,N	592	27	71	[6, 60]	[−3.62, 6.84]
LP	C,O,N,S,Cl	779	32	78	[6, 74]	[−3.62, 6.84]
SL	C,O,N	640	25	111	[4, 55]	[−9.33, 1.11]
SL	C,O,N,S,Cl	847	31	144	[4, 55]	[−11.60, 1.11]

shows the size and range of data sets that we prepared for each chemical property in Stage 1, where we denote the following:

• $\mathsf{\Lambda }$: the set of selected chemical elements (hydrogen atoms are added at the final stage);
• $|D_{\pi }|$: the size of data set $D_{\pi }$ over $\mathsf{\Lambda }$ for property $\pi$;
• $|{\mathsf{\Gamma }}^{\mathrm{int}}\left(D_{\pi }\right)|$: the number of different edge-configurations of interior-edges over the compounds in $D_{\pi }$;
• $\left|\mathcal{F}\right(D_{\pi }\left)\right|$: the number of non-isomorphic chemical rooted trees in the set of all 2-fringe-trees in the compounds in $D_{\pi }$;
• $[\underline{n},\overline{n}]$: the minimum and maximum values of $n\left(G\right)$ over the compounds G in $D_{\pi }$; and
• $[\underline{a},\overline{a}]$: the minimum and maximum values of $a\left(G\right)$ in $\pi$ over compounds G in $D_{\pi }$.

Stage 2. We used the new feature function that consists of the descriptors such as fringe-configuration defined in Section 2.1 and let $f_{\mathrm{fc}}$ denote the feature function.

Stage 3. Let $\eta :{\mathbb{R}}^K\to \mathbb{R}$ be a prediction function to a property function $a:D\to \mathbb{R}$ with a feature function $f:D\to {\mathbb{R}}^K$ for a data set D of chemical graphs. We define the coefficient of determination ${\mathrm{R}}^2(f,\eta ,D)$ of a prediction function $\eta$ over a data set D to be

$${\displaystyle {\mathrm{R}}^2(f,\eta ,D)\triangleq 1-\frac{{\sum }_{G\in D}{(a\left(G\right)-\eta \left(f\left(G\right)\right))}^2}{{\sum }_{G\in D}{(a\left(G\right)-\tilde{a})}^2}for\tilde{a}=\frac{1}{\left|D\right|}\sum _{G\in D}a\left(G\right).}$$

To conduct an experiment in Stage 3, we first constructed ten architectures $A_j$, $j\in [1,10]$ with one or two hidden layers. For each pair $(\pi ,A_j)$ of a property $\pi \in \{$K_OW, B_P, M_P, F_P, L_P, S_L}, and an architecture $A_j$, $j\in [1,10]$, we constructed five prediction functions in order to evaluate the performance with cross-validation as follows. Partition data set $D_{\pi }$ into five subsets $D_{\pi }^{\left(i\right)}$, $i\in [1,5]$ randomly and for each set $D_{\pi }\backslash D_{\pi }^{\left(i\right)}$ construct an ANN $\mathcal{N}(j,i)$ and its prediction function ${\eta }_{\mathcal{N}(j,i)}$ using the feature function $f_{\mathrm{fc}}$. We used scikit-learn version 0.23.2 with Python 3.8.5, MLPRegressor and ReLU activation function to construct each ANN $\mathcal{N}(j,i)$. We evaluated the resulting prediction function ${\eta }_{\mathcal{N}(j,i)}$ with the coefficient ${\mathrm{R}}^2(f_{\mathrm{fc}},{\eta }_{\mathcal{N}(j,i)},D_{\pi }^{\left(i\right)})$ of determination for the test set $D_{\pi }^{\left(i\right)}$. For each property $\pi$, let t-${\mathrm{R}}_{\mathrm{cv}}^2\left(j\right)$ denote the average of ${\mathrm{R}}^2(f_{\mathrm{fc}},{\eta }_{\mathcal{N}(j,i)},D_{\pi }^{\left(i\right)})$ over all $i\in [1,5]$ in the cross-validation to an architecture $A_j$.

Table 4. Results of Stages 2 and 3 in Phase 1.

$\mathit{\pi }$	$\mathsf{\Lambda }$	L-Time	t-${\mathbf{R}}_{\mathbf{cv}}^2$ (Best)	t-${\mathbf{R}}_{\mathbf{max}}^2$	Arch.
KOW	C,O,N	0.7	0.959	0.983	(156,10,10,1)
KOW	C,O,N,S,Cl	0.7	0.947	0.968	(199,20,10,1)
BP	C,O,N	3.5	0.858	0.923	(135,30,20,1)
BP	C,O,N,S,Cl	3.3	0.821	0.899	(163,10,1)
MP	C,O,N	3.8	0.784	0.893	(139,40,1)
MP	C,O,N,S,Cl	4.1	0.796	0.880	(170,10,10,1)
FP	C,O,N	1.1	0.750	0.874	(128,40,1)
FP	C,O,N,S,Cl	1.8	0.707	0.853	(157,10,10,1)
LP	C,O,N	0.5	0.868	0.908	(121,30,1)
LP	C,O,N,S,Cl	0.7	0.861	0.892	(137,20,10,1)
SL	C,O,N	0.7	0.870	0.913	(159,30,1)
SL	C,O,N,S,Cl	0.9	0.870	0.903	(201,30,20,1)

shows the results on Stages 2 and 3, where we denote the following.

-

$\mathsf{\Lambda }$: the set of selected chemical elements (hydrogen atoms are added at the final stage);

-

L-time: the average time (s) to construct an ANN over all $10\times 5=50$ ANNs;

-

t-${\mathrm{R}}_{\mathrm{cv}}^2$ (best): the best value of t-${\mathrm{R}}_{\mathrm{cv}}^2\left(j\right)$ over all architectures $A_j$, $j\in [1,10]$;

-

t-${\mathrm{R}}_{\max }^2$: the maximum of ${\mathrm{R}}^2(f_{\mathrm{fc}},{\eta }_{\mathcal{N}(j,i)},D_{\pi }^{\left(i\right)})$ over all $j\in [1,10],i\in [1,5]$; and

-

Arch.: The architecture $A_j$, $j\in [1,10]$ that attains t-${\mathrm{R}}_{\max }^2$. An architecture $(K,p,1)$ (resp., $(K,p_1,p_2,1)$) consists of an input layer with K nodes, a hidden layer with p nodes (resp., two hidden layers with $p_1$ and $p_2$ nodes, respectively), and an output layer with a single node, where K is equal to the number of descriptors in the feature vector.

From Table 4, we see that the execution of Stage 3 was considerably successful, where most of t-${\mathrm{R}}_{\max }^2$ are around 0.85 to 0.95 for all six chemical properties.

An Additional Experiment in Stage 3. We conducted an additional experiment to compare our new feature function $f_{\mathrm{fc}}$ with the feature function $f_{\mathrm{ec}}$ based edge-configuration in the previous method [27] designed with the same framework. Note that the previous feature vector $f_{\mathrm{ec}}\left(G\right)$ can be defined only for a cyclic graph G, whereas our feature vector $f_{\mathrm{fc}}\left(G\right)$ is defined for an arbitrary graph G. For each property $\pi \in \{$K_OW, B_P, M_P, F_P, L_P, S_L}, we set a set $\mathsf{\Lambda }$ of chemical elements to be $\left\{\mathrm{C},\mathrm{O},\mathrm{N},\mathrm{S},\mathrm{Cl}\right\}$ and then collect a data set ${\tilde{D}}_{\pi }$ of chemical cyclic graphs from the data set $D_{\pi }$ of all chemical graphs over the set $\mathsf{\Lambda }$ of chemical elements in the previous experiment. For each of the feature functions $f_{\mathrm{ec}}$ and $f_{\mathrm{fc}}$, we constructed five prediction functions with the same set of ten architectures $A_j$, $j\in [1,10]$ and the data set ${\tilde{D}}_{\pi }$ of chemical cyclic graphs in the same manner of the previous experiment.

Table 5. Results of prediction functions by $f_{\mathrm{ec}}$ and $f_{\mathrm{fc}}$ in data set ${\tilde{D}}_{\pi }$ of cyclic graphs and $f_{\mathrm{fc}}$ in data set $D_{\pi }$ of all graphs.

	$\mathit{f}={\mathit{f}}_{\mathbf{ec}}$, $\mathit{D}={\tilde{\mathit{D}}}_{\mathit{\pi }}$			$\mathit{f}={\mathit{f}}_{\mathbf{fc}}$, $\mathit{D}={\tilde{\mathit{D}}}_{\mathit{\pi }}$		$\mathit{f}={\mathit{f}}_{\mathbf{fc}}$, $\mathit{D}={\mathit{D}}_{\mathit{\pi }}$
$\mathit{\pi }$	$|{\tilde{\mathit{D}}}_{\mathit{\pi }}|$	t-${\mathbf{R}}_{\mathbf{cv}}^2$ (ave.)	t-${\mathbf{R}}_{\mathbf{cv}}^2$ (Best)	t-${\mathbf{R}}_{\mathbf{cv}}^2$ (ave.)	t-${\mathbf{R}}_{\mathbf{cv}}^2$ (Best)	$|{\mathit{D}}_{\mathit{\pi }}|$	t-${\mathbf{R}}_{\mathbf{cv}}^2$ (ave.)	t-${\mathbf{R}}_{\mathbf{cv}}^2$ (Best)
KOW	580	0.952	0.959	0.950	0.954	837	0.944	0.947
BP	224	0.688	0.718	0.680	0.693	425	0.809	0.821
MP	348	0.668	0.694	0.712	0.736	548	0.776	0.796
FP	218	0.435	0.476	0.574	0.623	399	0.688	0.707
LP	776	0.832	0.842	0.853	0.861	779	0.854	0.861
SL	638	0.851	0.863	0.853	0.861	847	0.860	0.870

shows the results of this experiment, where the table also includes the result of prediction functions by $f_{\mathrm{fc}}$ in the set $D_{\pi }$ of all chemical graphs. In the table, we denote the following:

-

$|{\tilde{D}}_{\pi }|$, $|D_{\pi }|$: the size of data set ${\tilde{D}}_{\pi }$ of cyclic graphs (resp., $D_{\pi }$ of all chemical graphs) for property $\pi$;

-

t-${\mathrm{R}}_{\mathrm{cv}}^2$ (ave.): the average of ${\mathrm{R}}^2(f,{\eta }_{\mathcal{N}(j,i)},D^{\left(i\right)})$ over all $j\in [1,10],i\in [1,5]$ for $f=f_{\mathrm{ec}},f_{\mathrm{fc}}$ and $D={\tilde{D}}_{\pi },D_{\pi }$; and

-

t-${\mathrm{R}}_{\mathrm{cv}}^2$ (best): ${max}_{j\in [1,10]}\{$the average of ${\mathrm{R}}^2(f_{\mathrm{fc}},{\eta }_{\mathcal{N}(j,i)},D_{\pi }^{\left(i\right)})$ over all $i\in [1,5]\}$.

From Table 5, we see that the score of R${}^2$ of the prediction function by $f_{\mathrm{fc}}$ in chemical cyclic graphs (resp., in all chemical graphs) is improved from that by $f_{\mathrm{ec}}$ for properties M_P and F_P (resp., B_P, M_P, and F_P). Recall that our new feature function $f_{\mathrm{fc}}$ can be defined for arbitrary graphs and we can select a larger data set than that by $f_{\mathrm{ec}}$ in a learning stage. This advantage is observed in the experiment. We guess that the better prediction function for B_P (resp., F_P) is obtained by using $f_{\mathrm{fc}}$ because the size of data set becomes considerably larger from $|{\tilde{D}}_{\pi }|=224$ to $|D_{\pi }|=425$ (resp., from $|{\tilde{D}}_{\pi }|=218$ to $|D_{\pi }|=399$).

Results on Phase 2.

We prepared the following instances (a–d) for conducting experiments of Stages 4 and 5 in Phase 2.

(a)

$I_{\mathrm{a}}=(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$: The instance used in Section 2.2 to explain the target specification.

(b)

$I_{\mathrm{b},i}=(G_{\mathrm{C}}^i,{\sigma }_{\mathrm{int}}^i,{\sigma }_{\mathrm{ce}}^i)$, $i=1,2,3,4$: An instance for inferring chemical graphs with rank at most 2. In the four instances $I_{\mathrm{b},i}$, $i=1,2,3,4$, the following specifications in $({\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ are common.

• Set $\mathsf{\Lambda }:=\left\{\mathrm{C},\mathrm{N},\mathrm{O}\right\}$, set ${\mathsf{\Lambda }}_{\mathrm{dg}}^{\mathrm{int}}$ to be the set of all possible symbols in $\mathsf{\Lambda }\times [1,4]$, and set ${\mathsf{\Gamma }}^{\mathrm{int}}$ to be the set of all possible edge-configurations. Set ${\mathsf{\Lambda }}^{*}\left(v\right):=\mathsf{\Lambda }$, $v\in V_{\mathrm{C}}$.
• The lower bounds ${\ell }_{\mathrm{LB}}$, ${\mathrm{bl}}_{\mathrm{LB}}$, ${\mathrm{ch}}_{\mathrm{LB}}$, ${\mathrm{bd}}_{2,\mathrm{LB}}$, ${\mathrm{bd}}_{3,\mathrm{LB}}$, ${\mathrm{na}}_{\mathrm{LB}}$, ${\mathrm{na}}_{\mathrm{LB}}^{\mathrm{int}}$, ${\mathrm{ns}}_{\mathrm{LB}}^{\mathrm{int}}$, ${\mathrm{ac}}_{\mathrm{LB}}^{\mathrm{int}}$, ${\mathrm{ec}}_{\mathrm{LB}}^{\mathrm{int}}$ are all set to be 0.
• The upper bounds ${\ell }_{\mathrm{UB}}$, ${\mathrm{bl}}_{\mathrm{UB}}$, ${\mathrm{ch}}_{\mathrm{UB}}$, ${\mathrm{bd}}_{2,\mathrm{UB}}$, ${\mathrm{bd}}_{3,\mathrm{UB}}$, ${\mathrm{na}}_{\mathrm{UB}}$, ${\mathrm{na}}_{\mathrm{UB}}^{\mathrm{int}}$, ${\mathrm{ns}}_{\mathrm{UB}}^{\mathrm{int}}$, ${\mathrm{ac}}_{\mathrm{UB}}^{\mathrm{int}}$, ${\mathrm{ec}}_{\mathrm{UB}}^{\mathrm{int}}$ are all set to be an upper bound $n^{*}$ on $n\left(G^{*}\right)$.
• For each property $\pi$, let $\mathcal{F}\left(D_{\pi }\right)$ denote the set of 2-fringe-trees in the compounds in $D_{\pi }$, and select a subset ${\mathcal{F}}_{\pi }^i\subseteq \mathcal{F}\left(D_{\pi }\right)$ with $|{\mathcal{F}}_{\pi }^i|=45-5i$, $i\in [1,5]$. For each instance $I_{\mathrm{b},i}$, set ${\mathcal{F}}_E:=\mathcal{F}\left(v\right):={\mathcal{F}}_{\pi }^i$, $v\in V_{\mathrm{C}}$.

Instance $I_{\mathrm{b},1}$ is given by the rank-1 seed graph $G_{\mathrm{C}}^1$ in Figure 9a and Instances $I_{\mathrm{b},i}$, $i=2,3,4$ are given by the rank-2 seed graph $G_{\mathrm{C}}^i$, $i=2,3,4$ in Figure 9b–d.

(i)

For instance $I_{\mathrm{b},1}$, select as a seed graph the monocyclic graph $G_{\mathrm{C}}^1=(V_{\mathrm{C}},E_{\mathrm{C}}=E_{(\ge 2)}\cup E_{(\ge 1)})$ in Figure 9a, where $V_{\mathrm{C}}=\{u_1,u_2\}$, $E_{(\ge 2)}=\left\{a_1\right\}$ and $E_{(\ge 1)}=\left\{a_2\right\}$. Set ${\mathrm{n}}_{\mathrm{LB}}^{\mathrm{int}}:=0,{\mathrm{n}}_{\mathrm{UB}}^{\mathrm{int}}:=12$ and $n_{\mathrm{LB}}:=n^{*}:=38$. We include a linear constraint $\ell \left(a_1\right)\le \ell \left(a_2\right)$ as part of the side constraint.

(ii)

For instance $I_{\mathrm{b},2}$, select as a seed graph the graph $G_{\mathrm{C}}^2=(V_{\mathrm{C}},E_{\mathrm{C}}=E_{(\ge 2)}\cup E_{(\ge 1)}\cup E_{(=1)})$ in Figure 9b, where $V_{\mathrm{C}}=\{u_1,u_2,u_3,u_4\}$, $E_{(\ge 2)}=\{a_1,a_2\}$, $E_{(\ge 1)}=\left\{a_3\right\}$ and $E_{(=1)}=\{a_4,a_5\}$. Set ${\mathrm{n}}_{\mathrm{LB}}^{\mathrm{int}}:={\mathrm{n}}_{\mathrm{UB}}^{\mathrm{int}}:=30$ and $n_{\mathrm{LB}}:=n^{*}:=50$. We include a linear constraint $\ell \left(a_1\right)\le \ell \left(a_2\right)$.

(iii)

For instance $I_{\mathrm{b},3}$, select as a seed graph the graph $G_{\mathrm{C}}^3=(V_{\mathrm{C}},E_{\mathrm{C}}=E_{(\ge 2)}\cup E_{(\ge 1)}\cup E_{(=1)})$ in Figure 9c, where $V_{\mathrm{C}}=\{u_1,u_2,u_3,u_4\}$, $E_{(\ge 2)}=\left\{a_1\right\}$, $E_{(\ge 1)}=\{a_2,a_3\}$ and $E_{(=1)}=\{a_4,a_5\}$. Set ${\mathrm{n}}_{\mathrm{LB}}^{\mathrm{int}}:={\mathrm{n}}_{\mathrm{UB}}^{\mathrm{int}}:=30$ and $n_{\mathrm{LB}}:=n^{*}:=50$. We include linear constraints $\ell \left(a_1\right)\le \ell \left(a_2\right)+\ell \left(a_3\right)$ and $\ell \left(a_2\right)\le \ell \left(a_3\right)$.

(iv)

For instance $I_{\mathrm{b},4}$, select as a seed graph the graph $G_{\mathrm{C}}^4=(V_{\mathrm{C}},E_{\mathrm{C}}=E_{(\ge 2)}\cup E_{(\ge 1)}\cup E_{(=1)})$ in Figure 9d, where $V_{\mathrm{C}}=\{u_1,u_2,u_3,u_4\}$, $E_{(\ge 1)}=\{a_1,a_2,a_3\}$ and $E_{(=1)}=\{a_4,a_5\}$. Set ${\mathrm{n}}_{\mathrm{LB}}^{\mathrm{int}}:={\mathrm{n}}_{\mathrm{UB}}^{\mathrm{int}}:=30$ and $n_{\mathrm{LB}}:=n^{*}:=50$. We include linear constraints $\ell \left(a_2\right)\le \ell \left(a_1\right)+1$, $\ell \left(a_2\right)\le \ell \left(a_3\right)+1$ and $\ell \left(a_1\right)\le \ell \left(a_3\right)$.

We define instances in (c) and (d) in order to find chemical graphs that have an intermediate structure of given two chemical cyclic graphs $G_A=(H_A=(V_A,E_A),{\alpha }_A,{\beta }_A)$ and $G_B=(H_B=(V_B,E_B),{\alpha }_B,{\beta }_B)$. Let ${\mathsf{\Lambda }}_A^{\mathrm{int}}$ and ${\mathsf{\Lambda }}_{\mathrm{dg},A}^{\mathrm{int}}$ denote the sets of chemical elements and chemical symbols of the interior-vertices in $G_A$, ${\mathsf{\Gamma }}_A^{\mathrm{int}}$ denote the sets of edge-configurations of the interior-edges in $G_A$, and ${\mathcal{F}}_A$ denote the set of 2-fringe-trees in $G_A$. Analogously define sets ${\mathsf{\Lambda }}_B^{\mathrm{int}}$, ${\mathsf{\Lambda }}_{\mathrm{dg},B}^{\mathrm{int}}$, ${\mathsf{\Gamma }}_B^{\mathrm{int}}$, and ${\mathcal{F}}_B$ in $G_B$.

(c)

$I_{\mathrm{c}}=(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$: An instance aimed to infer a chemical graph $G^{†}$ such that the core of $G^{†}$ is equal to the core of $G_A$ and the frequency of each edge-configuration in the non-core of $G^{†}$ is equal to that of $G_B$. We use chemical compounds CID 24822711 and CID 59170444 in Figure 10a,b for $G_A$ and $G_B$, respectively.

Set a seed graph $G_{\mathrm{C}}=(V_{\mathrm{C}},E_{\mathrm{C}}=E_{(=1)})$ to be the core of $G_A$.

Set $\mathsf{\Lambda }:=\left\{\mathrm{C},\mathrm{N},\mathrm{O}\right\}$, and set ${\mathsf{\Lambda }}_{\mathrm{dg}}^{\mathrm{int}}$ to be the set of all possible chemical symbols in $\mathsf{\Lambda }\times [1,4]$.

Set ${\mathsf{\Gamma }}^{\mathrm{int}}:={\mathsf{\Gamma }}_A^{\mathrm{int}}\cup {\mathsf{\Gamma }}_B^{\mathrm{int}}$ and ${\mathsf{\Lambda }}^{*}\left(v\right):=\left\{{\alpha }_A\left(v\right)\right\}$, $v\in V_{\mathrm{C}}$.

Set ${\mathrm{n}}_{\mathrm{LB}}^{\mathrm{int}}:=min\{{\mathrm{n}}^{\mathrm{int}}\left(G_A\right),{\mathrm{n}}^{\mathrm{int}}\left(G_B\right)\}$, ${\mathrm{n}}_{\mathrm{UB}}^{\mathrm{int}}:=max\{{\mathrm{n}}^{\mathrm{int}}\left(G_A\right),{\mathrm{n}}^{\mathrm{int}}\left(G_B\right)\}$,

$n_{\mathrm{LB}}:=min\{n\left(G_A\right),n\left(G_B\right)\}-10$ and $n^{*}:=max\{n\left(G_A\right),n\left(G_B\right)\}+5$.

Set lower bounds ${\ell }_{\mathrm{LB}}$, ${\mathrm{bl}}_{\mathrm{LB}}$, ${\mathrm{ch}}_{\mathrm{LB}}$, ${\mathrm{bd}}_{2,\mathrm{LB}}$, ${\mathrm{bd}}_{3,\mathrm{LB}}$, ${\mathrm{na}}_{\mathrm{LB}}$, ${\mathrm{na}}_{\mathrm{LB}}^{\mathrm{int}}$, ${\mathrm{ns}}_{\mathrm{LB}}^{\mathrm{int}}$ and ${\mathrm{ac}}_{\mathrm{LB}}^{\mathrm{int}}$ to be 0.

Set upper bounds ${\ell }_{\mathrm{UB}}$, ${\mathrm{bl}}_{\mathrm{UB}}$, ${\mathrm{ch}}_{\mathrm{UB}}$, ${\mathrm{bd}}_{2,\mathrm{UB}}$, ${\mathrm{bd}}_{3,\mathrm{UB}}$, ${\mathrm{na}}_{\mathrm{UB}}$, ${\mathrm{na}}_{\mathrm{UB}}^{\mathrm{int}}$, ${\mathrm{ns}}_{\mathrm{UB}}^{\mathrm{int}}$ and ${\mathrm{ac}}_{\mathrm{UB}}^{\mathrm{int}}$ to be $n^{*}$.

Set ${\mathrm{ec}}_{\mathrm{LB}}^{\mathrm{int}}\left(\gamma \right)$ to be the number of core-edges in $G_A$ with $\gamma \in {\mathsf{\Gamma }}^{\mathrm{int}}$ and ${\mathrm{ec}}_{\mathrm{UB}}^{\mathrm{int}}\left(\gamma \right)$ to be the number interior-edges in $G_A$ and $G_B$ with edge-configuration $\gamma$.

Let ${\mathcal{F}}_B^{\left(p\right)},p\in [1,2]$ denote the set of chemical rooted trees r-isomorphic p-fringe-trees in $G_B$.

Set ${\mathcal{F}}_E:=\mathcal{F}\left(v\right):={\mathcal{F}}_B^{\left(1\right)}\cup {\mathcal{F}}_B^{\left(2\right)}$, $v\in V_{\mathrm{C}}$.

(d)

$I_{\mathrm{d}}=(G_{\mathrm{C}}^1,{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$: An instance aimed to infer a chemical monocyclic graph $G^{†}$ such that the frequency vector of edge-configurations in $G^{†}$ is a vector obtained by merging those of $G_A$ and $G_B$. We use chemical monocyclic compounds CID 10076784 and CID 44340250 in Figure 10c,d for $G_A$ and $G_B$, respectively. Set a seed graph to be the monocyclic seed graph $G_{\mathrm{C}}^1=(V_{\mathrm{C}},E_{\mathrm{C}}=E_{(\ge 2)}\cup E_{(\ge 1)})$ with $V_{\mathrm{C}}=\{u_1,u_2\}$, $E_{(\ge 2)}=\left\{a_1\right\}$ and $E_{(\ge 1)}=\left\{a_2\right\}$ in Figure 9a.

Set $\mathsf{\Lambda }:=\left\{\mathrm{C},\mathrm{N},\mathrm{O}\right\}$, ${\mathsf{\Lambda }}_{\mathrm{dg}}^{\mathrm{int}}:={\mathsf{\Lambda }}_{\mathrm{dg},A}^{\mathrm{int}}\cup {\mathsf{\Lambda }}_{\mathrm{dg},B}^{\mathrm{int}}$ and ${\mathsf{\Gamma }}^{\mathrm{int}}:={\mathsf{\Gamma }}_A^{\mathrm{int}}\cup {\mathsf{\Gamma }}_B^{\mathrm{int}}$.

Set ${\mathrm{n}}_{\mathrm{LB}}^{\mathrm{int}}:=min\{{\mathrm{n}}^{\mathrm{int}}\left(G_A\right),{\mathrm{n}}^{\mathrm{int}}\left(G_B\right)\}$, ${\mathrm{n}}_{\mathrm{UB}}^{\mathrm{int}}:=max\{{\mathrm{n}}^{\mathrm{int}}\left(G_A\right),{\mathrm{n}}^{\mathrm{int}}\left(G_B\right)\}$,

$n_{\mathrm{LB}}:=min\{n\left(G_A\right),n\left(G_B\right)\}$ and $n^{*}:=max\{n\left(G_A\right),n\left(G_B\right)\}$.

Set lower bounds ${\ell }_{\mathrm{LB}}$, ${\mathrm{bl}}_{\mathrm{LB}}$, ${\mathrm{ch}}_{\mathrm{LB}}$, ${\mathrm{bd}}_{2,\mathrm{LB}}$, ${\mathrm{bd}}_{3,\mathrm{LB}}$, ${\mathrm{na}}_{\mathrm{LB}}$, ${\mathrm{na}}_{\mathrm{LB}}^{\mathrm{int}}$, ${\mathrm{ns}}_{\mathrm{LB}}^{\mathrm{int}}$ and ${\mathrm{ac}}_{\mathrm{LB}}^{\mathrm{int}}$ to be 0.

Set upper bounds ${\ell }_{\mathrm{UB}}$, ${\mathrm{bl}}_{\mathrm{UB}}$, ${\mathrm{ch}}_{\mathrm{UB}}$, ${\mathrm{bd}}_{2,\mathrm{UB}}$, ${\mathrm{bd}}_{3,\mathrm{UB}}$, ${\mathrm{na}}_{\mathrm{UB}}$, ${\mathrm{na}}_{\mathrm{UB}}^{\mathrm{int}}$, ${\mathrm{ns}}_{\mathrm{UB}}^{\mathrm{int}}$ and ${\mathrm{ac}}_{\mathrm{UB}}^{\mathrm{int}}$ to be $n^{*}$.

For each edge-configuration $\gamma \in {\mathsf{\Gamma }}^{\mathrm{int}}$, let ${\mathit{x}}_A^{*}\left({\gamma }^{\mathrm{int}}\right)$ (resp., ${\mathit{x}}_B^{*}\left({\gamma }^{\mathrm{int}}\right)$) denote the number of interior-edges with $\gamma$ in $G_A$ (resp., $G_B$), $\gamma \in {\mathsf{\Gamma }}^{\mathrm{int}}$ and set

${\mathit{x}}_{min}^{*}\left(\gamma \right):=min\{{\mathit{x}}_A^{*}\left(\gamma \right),{\mathit{x}}_B^{*}\left(\gamma \right)\}$, ${\mathit{x}}_{max}^{*}\left(\gamma \right):=max\{{\mathit{x}}_A^{*}\left(\gamma \right),{\mathit{x}}_B^{*}\left(\gamma \right)\}$,

${\mathrm{ec}}_{\mathrm{LB}}^{\mathrm{int}}\left(\gamma \right):=\lfloor (3/4){\mathit{x}}_{min}^{*}\left(\gamma \right)+(1/4){\mathit{x}}_{max}^{*}\left(\gamma \right)\rfloor$ and

${\mathrm{ec}}_{\mathrm{UB}}^{\mathrm{int}}\left(\gamma \right):=\lceil (1/4){\mathit{x}}_{min}^{*}\left(\gamma \right)+(3/4){\mathit{x}}_{max}^{*}\left(\gamma \right)\rceil$.

Set ${\mathcal{F}}_E:=\mathcal{F}\left(v\right):={\mathcal{F}}_A\cup {\mathcal{F}}_B$, $v\in V_{\mathrm{C}}$.

In Stage 5, before we formulate an MILP for inferring a target chemical graph $G^{†}$ for each instance I, we reduce the input layer of an ANN $\mathcal{N}$ constructed in Stage 3 so that the input layer consists of input nodes that correspond to the descriptors actually used in the specification $(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ of the instance I, i.e., we remove any input nodes in $\mathcal{N}$ that represent the frequency of edge-configurations in ${\mathsf{\Gamma }}^{\mathrm{int}}\left(D_{\pi }\right)$ and chemical rooted trees $\psi \in \mathcal{F}\left(D_{\pi }\right)$ not contained in the specification $(G_{\mathrm{C}},{\sigma }_{\mathrm{int}},{\sigma }_{\mathrm{ce}})$ of I. For example, there are $\left|\mathcal{F}\right(D_{\pi }\left)\right|=109$ chemical rooted trees in the set of 2-fringe-trees in the data set $D_{\pi }$ with $\pi =$ K_OW in Table 3, and an ANN $\mathcal{N}$ constructed in Stage 3 contains 109 input nodes that correspond to the descriptors for the fringe-configuration. However, the set of input nodes for the fringe-configuration is reduced to a set of $|{\mathcal{F}}^{*}|=40$ input nodes when we formulate an MILP for solving instance $I_{\mathrm{b},1}$, saving the number of integer variables.

Table 6. Features of test instances.

Instance	$\mathsf{\Lambda }$	$|{\mathsf{\Gamma }}^{\mathbf{int}}|$	$|{\mathcal{F}}^{*}|$	$[{\mathbf{n}}_{\mathbf{LB}}^{\mathbf{int}},{\mathbf{n}}_{\mathbf{UB}}^{\mathbf{int}}]$	$[{\mathit{n}}_{\mathbf{LB}},{\mathit{n}}^{*}]$
$I_{\mathrm{a}}$	C,O,N	10	11	[30,50]	[20,28]
$I_{\mathrm{b},1}$	C,O,N	28	40	[38,38]	[6,6]
$I_{\mathrm{b},2}$	C,O,N	28	35	[50,50]	[30,30]
$I_{\mathrm{b},3}$	C,O,N	28	30	[50,50]	[30,30]
$I_{\mathrm{b},4}$	C,O,N	28	25	[50,50]	[30,30]
$I_{\mathrm{c}}$	C,O,N	8	12	[46,46]	[24,24]
$I_{\mathrm{d}}$	C,O,N	7	8	[40,45]	[18,18]

shows the features of the seven test instances, where we denote the following:

-

$\mathsf{\Lambda }$: the set of non-hydrogen chemical elements for inferring a target graph;

-

$|{\mathsf{\Gamma }}^{\mathrm{int}}|$: the number of different edge-configurations of interior-edges for inferring a target graph;

-

$|{\mathcal{F}}^{*}|$: the number of different chemical rooted trees in the set ${\mathcal{F}}^{*}={\mathcal{F}}_E\cup {\bigcup }_{v\in V_{\mathrm{C}}}\mathcal{F}\left(v\right)$; and

-

$[{\mathrm{n}}_{\mathrm{LB}}^{\mathrm{int}},{\mathrm{n}}_{\mathrm{UB}}^{\mathrm{int}}]$, $[n_{\mathrm{LB}},n^{*}]$: the lower and upper bounds on ${\mathrm{n}}^{\mathrm{int}}\left(G^{†}\right)$ and $n\left(G^{†}\right)$ for inferring a target graph $G^{†}$.

Stage 4. To solve an MILP in Stage 4, we used CPLEX version 12.10.

Table 7. Results of Stage 4 for K_OW.

Instance	$[\underline{\mathit{a}},\overline{\mathit{a}}]$	$[\underline{\mathit{y}},\overline{\mathit{y}}]$	${\mathit{y}}^{*}$	#v	#c	IP-Time	n	${\mathbf{n}}^{\mathbf{int}}$
$I_{\mathrm{a}}$	[−7.53, 13.45]	[−7.0, 13.4]	3.2	7663	9162	3.9	35	24
$I_{\mathrm{b},1}$	[−7.53, 13.45]	[−7.5, 13.4]	3.0	9894	6626	17.5	38	7
$I_{\mathrm{b},2}$	[−7.53, 13.45]	[−7.5, 13.4]	3.0	11,514	8934	14.0	50	30
$I_{\mathrm{b},3}$	[−7.53, 13.45]	[−7.5, 13.4]	3.0	11,318	8926	24.6	50	30
$I_{\mathrm{b},4}$	[−7.53, 13.45]	[−7.5, 13.4]	3.0	11,122	8918	22.0	50	30
$I_{\mathrm{c}}$	[−7.53, 13.45]	[−7.5, 13.4]	3.0	7867	8630	2.1	49	32
$I_{\mathrm{d}}$	[−7.53, 13.45]	[−7.5, 13.4]	3.0	5395	6899	5.2	45	23

,

Table 8. Results of Stage 4 for B_P.

Instance	$[\underline{\mathit{a}},\overline{\mathit{a}}]$	$[\underline{\mathit{y}},\overline{\mathit{y}}]$	${\mathit{y}}^{*}$	#v	#c	IP-Time	n	${\mathbf{n}}^{\mathbf{int}}$
$I_{\mathrm{a}}$	[−11.70, 470.0]	[352, 470]	411	7583	8982	2.7	42	25
$I_{\mathrm{b},1}$	[−11.70, 470.0]	[−11, 470]	229	9816	6449	2.7	38	7
$I_{\mathrm{b},2}$	[−11.70, 470.0]	[−11, 470]	229	11,436	8757	9.1	50	30
$I_{\mathrm{b},3}$	[−11.70, 470.0]	[−11, 470]	229	11,240	8749	11.0	50	30
$I_{\mathrm{b},4}$	[−11.70, 470.0]	[−11, 470]	229	11,044	8741	24.0	50	30
$I_{\mathrm{c}}$	[−11.70, 470.0]	[170, 470]	320	7575	8450	25.9	49	33
$I_{\mathrm{d}}$	[−11.70, 470.0]	[151, 470]	310	5315	6719	4.4	43	23

,

Table 9. Results of Stage 4 for M_P.

Instance	$[\underline{\mathit{a}},\overline{\mathit{a}}]$	$[\underline{\mathit{y}},\overline{\mathit{y}}]$	${\mathit{y}}^{*}$	#v	#c	IP-Time	n	${\mathbf{n}}^{\mathbf{int}}$
$I_{\mathrm{a}}$	[−185.3, 300.0]	[55, 300]	177.5	7602	9023	16.1	41	24
$I_{\mathrm{b},1}$	[−185.3, 300.0]	[−180, 300]	60	9833	6487	2.3	38	9
$I_{\mathrm{b},2}$	[−185.3, 300.0]	[−185, 300]	57.4	11,453	8795	44.7	50	30
$I_{\mathrm{b},3}$	[−185.3, 300.0]	[−185, 300]	57.4	11,257	8787	10.5	50	30
$I_{\mathrm{b},4}$	[−185.3, 300.0]	[−185, 300]	57.4	11,061	8779	93.9	50	30
$I_{\mathrm{c}}$	[−185.3, 300.0]	[253, 300]	260.0	7580	6172	24.0	41	33
$I_{\mathrm{d}}$	[−185.3, 300.0]	[−75, 299]	58	5110	4050	104.6	45	23

,

Table 10. Results of Stage 4 for F_P.

Instance	$[\underline{\mathit{a}},\overline{\mathit{a}}]$	$[\underline{\mathit{y}},\overline{\mathit{y}}]$	${\mathit{y}}^{*}$	#v	#c	IP-Time	n	${\mathbf{n}}^{\mathbf{int}}$
$I_{\mathrm{a}}$	[−82.99, 300.0]	[98, 300]	199	7459	8696	1.6	35	22
$I_{\mathrm{b},1}$	[−82.99, 300.0]	[−82, 300]	109	9694	6166	1.4	38	8
$I_{\mathrm{b},2}$	[−82.99, 300.0]	[−82, 300]	109	11,314	8474	8.7	50	30
$I_{\mathrm{b},3}$	[−82.99, 300.0]	[−82, 300]	109	11,118	8466	25.8	50	30
$I_{\mathrm{b},4}$	[−82.99, 300.0]	[−82, 300]	109	10,922	8458	8.5	50	30
$I_{\mathrm{c}}$	[−82.99, 300.0]	[250, 300]	275	7667	8170	60.9	47	34
$I_{\mathrm{d}}$	[−82.99, 300.0]	[54, 300]	177	5193	6436	2.0	45	23

,

Table 11. Results of Stage 4 for L_P.

Instance	$[\underline{\mathit{a}},\overline{\mathit{a}}]$	$[\underline{\mathit{y}},\overline{\mathit{y}}]$	${\mathit{y}}^{*}$	#v	#c	IP-Time	n	${\mathbf{n}}^{\mathbf{int}}$
$I_{\mathrm{a}}$	[−3.6, 6.84]	[−3.6, 6.8]	1.6	7597	9008	1.9	39	23
$I_{\mathrm{b},1}$	[−3.6, 6.84]	[−3.6, 6.8]	1.6	9836	6481	2.9	38	8
$I_{\mathrm{b},2}$	[−3.6, 6.84]	[−3.6, 6.8]	1.6	11,456	8789	21.1	50	30
$I_{\mathrm{b},3}$	[−3.6, 6.84]	[−3.6, 6.8]	1.6	11,260	8781	20.4	50	30
$I_{\mathrm{b},4}$	[−3.6, 6.84]	[−3.6, 6.8]	1.6	11,064	8773	24.2	50	30
$I_{\mathrm{c}}$	[−3.6, 6.84]	[−3.6, 6.8]	1.6	7801	8476	1.1	47	32
$I_{\mathrm{d}}$	[−3.6, 6.84]	[−3.6, 6.8]	1.6	5335	6754	4.3	45	23

and

Table 12. Results of Stage 4 for S_L.

Instance	$[\underline{\mathit{a}},\overline{\mathit{a}}]$	$[\underline{\mathit{y}},\overline{\mathit{y}}]$	${\mathit{y}}^{*}$	#v	#c	IP-Time	n	${\mathbf{n}}^{\mathbf{int}}$
$I_{\mathrm{a}}$	[−9.33, 1.11]	[−9.3, −2.0]	−5.6	7674	9186	2.4	41	23
$I_{\mathrm{b},1}$	[−9.33, 1.11]	[−9.3, −2.0]	−5.6	9906	6650	22.3	38	12
$I_{\mathrm{b},2}$	[−9.33, 1.11]	[−9.3, −2.0]	−5.6	11,526	8958	15.2	50	30
$I_{\mathrm{b},3}$	[−9.33, 1.11]	[−9.3, −2.0]	−5.6	11,330	8950	16.2	50	30
$I_{\mathrm{b},4}$	[−9.33, 1.11]	[−9.3, −2.0]	−5.6	11,134	8942	122.7	50	30
$I_{\mathrm{c}}$	[−9.33, 1.11]	[−9.3, −2.0]	−5.6	7874	8648	1.2	54	33
$I_{\mathrm{d}}$	[−9.33, 1.11]	[−9.3, −3.0]	−6.1	5402	6917	8.1	43	23

show the results on Stages 4 and 5, where we denote the following:

-

$[\underline{a},\overline{a}]$: the minimum and maximum values of $a\left(G\right)$ in $\pi$ over compounds G in $D_{\pi }$ in Table 3;

-

$[\underline{y},\overline{y}]$: $\underline{y}$ (resp., $\overline{y}$) denotes the minimum (resp., maximum) target value y with $\underline{a}\le y\le \overline{a}$ such that the MILP instance for the target value $y^{*}=y$ becomes feasible (i.e., admits a target chemical graph $G^{†}$). To determine the minimum and minimum target values $\underline{y}$ and $\overline{y}$, we solved many numbers of MILP instances. Note that the MILP instance may become infeasible for some value y within the range $[\underline{y},\overline{y}]$;

-

$y^{*}$: a target value in $[\underline{y},\overline{y}]$ for a property $\pi$;

-

#v: the number of variables in the MILP in Stage 4;

-

#c: the number of constraints in the MILP in Stage 4;

-

IP-time: the time (sec.) to solve the MILP in Stage 4;

-

n: the number $n\left(G^{†}\right)$ of non-hydrogen atoms in the chemical graph $G^{†}$ inferred in Stage 4; and

-

${\mathrm{n}}^{\mathrm{int}}$: the number ${\mathrm{n}}^{\mathrm{int}}\left(G^{†}\right)$ of interior-vertices in the chemical graph $G^{†}$ inferred in Stage 4.

Figure 11a illustrates the chemical graph $G^{†}$ inferred from instance $I_{\mathrm{c}}$ with $y^{*}=3.0$ of K_OW in Table 7.

Figure 11b illustrates the chemical graph $G^{†}$ inferred from instance $I_{\mathrm{d}}$ with $y^{*}=1.6$ of L_P in Table 11.

The topological specification of instances $I_{\mathrm{a}}$, $I_{\mathrm{c}}$ and $I_{\mathrm{d}}$ is more restricted than that of the other instances, and thereby the feasible target range $[\underline{y},\overline{y}]$ of $I_{\mathrm{a}}$, $I_{\mathrm{c}}$ and $I_{\mathrm{d}}$ is rather narrower than the original range $[\underline{a},\overline{a}]$ for some property $\pi$. We see that the running time for solving an MILP instance with $n=50$ is 8.5 to 122 (s), which is much smaller than the running time of 61 to 12058 (s) to solve a similar set of MILP instances with $n=50$ in the experimental result for the previous method [28].

Stage 5. We computed chemical isomers $G^{*}$ of each target chemical graph $G^{†}$ inferred in Stage 4. We execute the algorithm for generating chemical isomers of $G^{†}$ up to 100 when the number of all chemical isomers exceeds 100. The algorithm can evaluate a lower bound on the total number of all chemical isomers $G^{†}$ without generating all of them.

Table 13. Results of Stage 5 for K_OW, L_P, and B_P.

		Kow			Lp			Bp
Instance	DP-Time	G-LB	#G	DP-Time	G-LB	#G	DP-time	G-LB	#G
$I_{\mathrm{a}}$	0.031	16	16	0.164	128	100	0.164	$1.4\times {10}^5$	100
$I_{\mathrm{b}}^1$	0.149	$2.8\times {10}^5$	100	0.148	$2.0\times {10}^{10}$	100	0.162	$4.4\times {10}^5$	100
$I_{\mathrm{b}}^2$	44.1	$3.9\times {10}^{10}$	100	118	900	100	171	6	6
$I_{\mathrm{b}}^3$	27.2	20	20	80.2	6	6	28.6	7	7
$I_{\mathrm{b}}^4$	0.166	6000	100	73	12	12	142	5	5
$I_{\mathrm{c}}$	0.166	6000	100	0.168	288	100	0.168	$4.0\times {10}^5$	100
$I_{\mathrm{d}}$	22.3	$8.3\times {10}^{10}$	100	1.44	$3.2\times {10}^8$	100	1.7	$9.7\times {10}^9$	100

and

Table 14. Results of Stage 5 for F_P, M_P, and S_L.

		FP			MP			SL
Instance	DP-Time	G-LB	#G	DP-Time	G-LB	#G	DP-Time	G-LB	#G
$I_{\mathrm{a}}$	0.057	32	32	0.165	256	100	0.165	1024	100
$I_{\mathrm{b}}^1$	0.164	$3.1\times {10}^6$	100	0.166	$1.4\times {10}^6$	100	0.163	$4.5\times {10}^5$	100
$I_{\mathrm{b}}^2$	28.8	720	100	8.26	$2.4\times {10}^{10}$	100	1.07	$5.6\times {10}^9$	100
$I_{\mathrm{b}}^3$	72.2	27	27	51.9	1	1	46.5	1680	100
$I_{\mathrm{b}}^4$	40.3	20	20	125	$6.1\times {10}^7$	100	7.01	$1.1\times {10}^8$	100
$I_{\mathrm{c}}$	0.169	$1.1\times {10}^5$	100	0.173	6048	100	0.168	120	100
$I_{\mathrm{d}}$	0.057	32	32	0.17	$4.2\times {10}^8$	100	0.165	1024	100

show the computational results of the experiment, where we denote the following:

-

DP-time: the running time (s) to execute the dynamic programming algorithm in Stage 5 to compute a lower bound on the number of all chemical isomers $G^{*}$ of $G^{†}$ and generate all (or up to 100) chemical isomers $G^{*}$;

-

G-LB: a lower bound on the number of all chemical isomers $G^{*}$ of $G^{†}$; and

-

#G: the number of all (or up to 100) chemical isomers $G^{*}$ of $G^{†}$ generated in Stage 5.

From Table 13 and Table 14, we observe that the running time for generating up to 100 target chemical graphs in Stage 5 is not considerably larger than that in Stage 4.

4. Discussions and Conclusions

The framework of designing chemical graphs using ANNs and MILP has been proposed [23] as a basis of a total system of the QSAR and the inverse of QSAR, where the inverse of a prediction function produced by an ANN is solved by an MILP. The merit of the framework is that the inverse problem can be treated exactly as a mathematical problem, and an MILP instance with a moderate size can be efficiently solved with a fast MILP solver. On the other hand, the main technical concern in applying the framework is in defining a feature vector of a chemical graph in terms of graph theoretical descriptors so that the computation of a feature vector can be simulated with a set of linear constraints in an MILP. So far, the framework has been applied to the design of new methods of inferring several restricted classes of chemical graphs such as the graphs with rank at most 2 and the $\rho$-lean cyclic graphs [26,28].

Herein, we examine some technical issues in the previous method before we observe some new features of our method in this paper.

In the feature vector of the previous models [26,28], the structure of subgraphs used as descriptors is only a pair of adjacent vertices, called adjacency-configuration or edge-configuration, which is significantly limited from a variety of subgraphs used in a more sophisticated construction of a feature vector such as the fingerprint. However, including the occurrence of a certain subgraph with only a few vertices as part of a feature vector may require realizing a mechanism of the subgraph isomorphism in an MILP that simulates the computation of such an occurrence and can easily make the resulting MILP very complicated and hard to solve. Furthermore, the feature vector can be defined only for cyclic graphs and we need to eliminate any acyclic graphs from the original data set before we construct a prediction function. This may reduce a data set to an unnecessarily small size or reduce the chances of capturing important information on QSAR over all types of graphs.

A branch-parameter $\rho$ was originally introduced as a new measure to the “agglomeration degree” of trees [24] and then used to define restricted classes of acyclic and cyclic graphs [24,27]. In fact, such a restriction on the structure of target chemical graphs was rather necessary to reduce the size of an MILP formulation that simulates a selection process of a target chemical graph from a supergraph (called a scheme graph), where the number of variables and constraints required to infer a chemical graph with $n^{*}$ non-hydrogen atoms is $O\left(n^{*}\right)$ when some other parameters such as $\rho$ are regarded as constants.

Although nearly 97% of cyclic chemical compounds with up to 100 non-hydrogen atoms in PubChem are 2-lean [24], the way of specifying the topological structure of a target chemical graph in the previous method [26,28] was based on the core and the non-core of a chemical graph, and we could not include a fixed substructure in the non-core of a target chemical graph.

Compared with the previous models, the two-layered model proposed in this paper is rather simple, where a chemical graph is regarded as a combination of the interior and the exterior. The new model can deal with chemical compounds with any graph structure and include a prescribed structure in both of the acyclic and cyclic parts of a target chemical graph as long as the requirement on target chemical graphs is described under the set of specification rules introduced in this paper. This considerably improves the availability of the framework in a practical application.

The feature vector of our two-layered model can be defined for arbitrary graphs. In the new feature vector, the exterior of a chemical graph is encoded into fringe-configurations, i.e., the occurrence of each chemical rooted tree with height at most $\rho$, where we may regard that the set of such a chemical rooted trees plays a similar role of some types of functional groups. In our method, we include as part of the descriptors of a feature vector the occurrence of each of such chemical rooted trees and the descriptors of our feature vector on the exterior of a chemical graph may have an analogous effect with the fingerprint.

Our specification of target chemical graphs can specify a candidate set $\mathcal{F}$ of chemical rooted trees that are allowed to be used as chemical rooted trees in the exterior of a target chemical graph. This allows us to control the chemical property of target chemical graphs in a more meaningful way since chemical properties of some rooted trees in $\mathcal{F}$ are known as functional groups and some kinds of rooted trees can be prohibited in a target chemical graph, if necessary, just by excluding from the candidate set $\mathcal{F}$. Although the number $\left|\mathcal{F}\right(D_{\pi }\left)\right|$ of different kinds of such chemical trees in a data set $D_{\pi }$ from PubChem is approximately up to 300 for $\rho =2$ in many cases and the number of input nodes in an ANN $\mathcal{N}$ becomes over $\left|\mathcal{F}\right(D_{\pi }\left)\right|$, we derived an MILP formulation for inferring a chemical graph with with $n^{*}$ non-hydrogen atoms and a candidate set $\mathcal{F}$ of chemical rooted trees by using $O(n^{*}+|\mathcal{F}\left|\right)$ variables and $O\left(n^{*}\right|\mathcal{F}\left|\right)$ constraints when the number of interior-vertices is constant, where $\left|\mathcal{F}\right|$ can be quite small compared with $\left|\mathcal{F}\right(D_{\pi }\left)\right|$.

We have implemented the proposed method for inferring chemical compounds with a prescribed topological substructure setting $\rho =2$. The results of computational experiments using some chemical properties such as octanol/water partition coefficient, boiling point, melting point, flash point, lipophylicity, and solubility suggest that the proposed system can infer chemical graphs with 50 non-hydrogen atoms.

For a larger branch-parameter, say $\rho =3,4$, we obtain a more variety of chemical rooted trees which provides new descriptors in a feature vector and new candidates for fringe-trees in the exterior in a target chemical graph, whereas the number of different chemical rooted trees in $\mathcal{F}\left(D_{\pi }\right)$ may increase rapidly.

It is left as a future work to use other learning methods such as decision tree, random forest, graph convolution, and an ensemble method in Stages 3 and 4 in the framework.